Quantifying economic risks to dairy farms from volcanic hazards in Taranaki, Aotearoa / New Zealand

The volcanic hazard and risk science for Taranaki Mounga (Taranaki volcano) in Aotearoa / New Zealand is in an advanced state, with robust probabilistic data and a series of direct impact scenarios modelled for the region. Here, we progress this work and demonstrate a method to provide risk information that is nuanced for factors such as location and economic sector and considers the dynamic nature of volcanism with hazards potentially repeated over time. Recognising the fundamental importance of the dairy sector to Taranaki region, this paper provides valuable insights into the potential impacts and risks to heterogeneous dairy cattle farms within the region from volcanic hazards. We provide volcanic impact and risk metrics in economic or monetary terms in order to improve its relevance to decision-makers while reducing the complexity of the impacts. To do this, we developed a dynamic, multi-event farm system model of response and recovery, which takes in hazard intensity metrics from a series of volcanic events and generates the resulting annualised revenues, expenditures, and recovery costs through time. The model is formulated in a generalised way such that it can be used for various other hazard types and agricultural land uses. In our application of the model, we create and apply a suite of 10 000 simulations that capture different iterations of possible future volcanic activity over a 50-year period. These include the generation of lahars following eruptions and associated failures for transport and water supply networks. Farms at five case study locations were modelled to capture the diversity in farm management and the spatial variation in hazard intensities and likelihoods across the region. We provide summaries of the distributions of economic impacts generated, both for individual events and for the 50-year volcanic future horizon. Drawing the information together, we also summarise the results for each case study farm in terms of the value at risk statistic. For the case study farms with negligible lahar risk, we find, with 90 % confidence, that volcanic losses over the next 50 years will not exceed around 10 % of property value. By comparison, for the farm with the most severe lahar and ashfall exposure, we find that, at the same level of confidence, losses extend to approximately half the property value. These results indicate that with access to sufficient risk information, we should anticipate volcanic risk as playing an important role in shaping the future dairy sector in Taranaki region. The modelling pipeline and assessment metrics demonstrated in this paper could be used to assess mitigation and adaptation strategies to reduce the risk from volcanic hazards and improve the resilience of farm businesses.

1 Introduction

Taranaki Mounga is a stratovolcano located in the western North Island of Aotearoa / New Zealand. The probability of a volcanic event occurring has been estimated at between 33 % and 42 % over the next 50 years (Damaschke et al., 2018). Thus, for the surrounding region, there is a distinct possibility of a volcanic event impacting economic activities, including the more than 1000 dairy farms in the vicinity of the mountain.

Although it is generally recognised that there is a smaller body of disaster risk science for volcanic hazards than other hazards, especially when compared to the more frequently occurring weather-related hazards (Ward et al., 2020), the hazard and risk assessment for Taranaki has a comparatively high level of research that has gone into understanding geological, engineering, and social aspects of risk. A body of research has helped to quantify eruption frequencies and magnitudes along with the variety of physical processes and consequential hazard types (see Cronin et al., 2021), while several studies have built knowledge of the vulnerabilities in socio-ecological systems (e.g. Wilson et al., 2009; Wild et al., 2019; Weir et al., 2024b) and started to quantify the impacts of potential eruption scenarios (Weir et al., 2022, 2024a; McDonald et al., 2017). Nevertheless, providing information on the dynamic impacts and risks to socio-economic activities that is nuanced for factors such as sector, location, timing, and management approach, as well as in a form that is appropriate to and meaningful for decision-makers, is an ongoing topic of research.

In this study, we focus on potential volcanic impacts and risks to dairy cattle farms within Taranaki region, recognising that dairy cattle farming is the dominant land use surrounding the volcano and that the dairy sector is strongly embedded in the regional economy directly constituting approximately 13 % of gross regional product (Cardwell et al., 2023). While it is feasible that over the short- to medium-term (say 20–50 years) no volcanic events will occur, it is also possible that there will be several events. We know that the avenues by which dairy farms may be impacted from volcanic events are various and complex, including multiple types of hazards (e.g. ashfall, lahars, volcanic mudflows, and floods) and hazards that manifest at the physical location of the farms as well as impacting on networks and infrastructure supporting farms, such as transport and water supply networks (Deligne et al., 2022). The effects from potentially multiple eruptive events or hazards may occur separately or overlap (Lawrence et al., 2020). Some impacts will be immediate and short term, while others will persist over time. Furthermore, heterogeneity across dairy farms (Doole and Pannell, 2012) and in human responses and adaptation will potentially influence the way in which impacts manifest.

Providing volcanic impact metrics in economic or monetary terms is one tactic to help reduce complexity and assist decision-makers in making sense of potential impact and risk information. Monetary metrics can aggregate across direct and indirect, immediate and future, and short- and long-term impacts. In this study, we propagate the impacts of volcanic events through individual farm system business models of response and recovery. The model evaluates economic impacts, considering both farm enterprise asset stocks and changes in expenditure and revenue flows. The use of such system-based models of enterprises is advantageous in terms of familiarity to stakeholders, having been used to quantify economic impacts for numerous policy and risk assessment contexts in Aotearoa / New Zealand, including the assessment of options to reduce agriculture greenhouse gas emissions and improve freshwater quality (see, for example, Parsons et al., 2015; Journeaux et al., 2020, and He Waka Eke Noa, 2022). Because the model provides a systemic view of farm operations, it is also possible to connect it to wider economic models that trace flow-on impacts through an economy (Monge et al., 2023; Hallegatte, 2008; Monge and McDonald, 2020; Brown et al., 2019a).

In summary, the purpose of this paper is to demonstrate how we can provide valuable insights into the potential economic impacts and risks to land-based enterprises around volcanic regions. To achieve this, we develop a dynamic business impact and recovery model that has broad applicability to volcanic as well as other hazard types, and we demonstrate its application in a modelling pipeline that quantifies economic impacts across many simulations and ultimately derives risk-based metrics. The example we have chosen as a case study is dairy cattle farms within Taranaki region in Aotearoa / New Zealand. The following sections of this paper begin by providing a general description of the dynamic business impact model developed for calculating economic impacts of hazard events at the farm scale (Sect. 2), as well as further details of the Taranaki case study (Sect. 3). This is followed (Sects. 4 and 5) by an explanation of the way in which the model was set up for the case study in terms of scenarios considered and data and parameter settings. In the “Results and discussion” section (Sect. 6), the economic impacts of individual volcanic events and over 50-year simulations are presented along with summaries of risk in terms of probability density functions and the value at risk statistic.

2 Agricultural Business Behaviours Model

The Agricultural Business Behaviours Model (AgriBBM) operates at a single farm level, inputting hazard intensity metrics (or metrics describing intensity of cascading failures such as water supply and transport disruption) from a series of events to work out the resulting annualised revenue, expenditure, and costs (REC) through time. In describing the model, we use notation in a way that we intend to be generalisable to multiple hazard types and agricultural land uses. The details developed for the specific case of dairy farms under volcanic hazards, including some necessary adjustments and extensions for the specific case study region, are covered in Sect. 4.3–4.6 and Appendix A to C.

Figure 1 Agricultural Business Behaviours Model overview diagram for a single farm, showing the progression from the hazard event information, HIM_e, and time, $t_e^{\mathrm{0}}$, for events e (orange) through to the results (green). Specifically, the annual revenue, expenditure, and costs, REC_i, for year i, which represents a sum of the adjusted farm accounts $\widehat{F_{i,j}}$, where j is the index for each line item of the annual revenue and expenditure, and the total annual recovery costs C_i. For each farm, the model needs the inputs (brown): base accounts at the line-item level for each year F_i,j for business as usual, base farm asset information (e.g. building replacement costs), and an estimate of the relative impact scalars at the line-item level, $S_j^{\mathrm{0}}$, for all possible overall combined impact states (ISs). All intermediary variables (blue and grey) are defined in the main text.

In this section, we will first step through the model for a single event, then describe the development of the agricultural business recovery trajectory function, and finally describe how we deal with multiple, potentially overlapping, events in a single scenario. An overview of the model is shown in Fig. 1.

2.1 AgriBBM steps for a single event

Each farm is modelled independently. For each farm, we start with a vector of base farm accounts F_i,j, where j is the index for each line item of the annual revenue and expenditure for business as usual, and i is the year. The annual farm operating profit is calculated from ${\sum }_j{F}_{i,j}$. In general, these accounts can be set to vary by year i. This allows for the inclusion of policies that change fees or subsidies through time or for farms where the revenue and expenditure depend on the age of the land use. This is relevant for agricultural industries such as horticulture where it can take a number of years for plants to become established, or forestry, which have very different annual expenditure and revenue depending on the age of the plantation.

As shown in Fig. 1, disruptions to the farm's operation and behaviour is included through scalars applied to base revenues, expenses, and recovery costs, in combination with a recovery trajectory determined through a recovery time period starting from the event at time t⁰.

First, we use the hazard intensity metrics to determine the damage states (DSs) for the event split up by each type of damage or disruption. The term “damage state” is chosen here simply for convenience. Note that we do allow for the individual categories and associated functions that could be considered here to include those relating to damages to capital as well as those describing the loss or disruption in the provision of a service, e.g. DS^ashfall for the damage state due to direct ashfall on the farm impacting on pastures, crops, stock, buildings, and other assets and DS^water for the damage state due to outages in the water supply to farms. Such damage states are typically determined by probabilistic fragility functions, and the set of damage states considered will depend on the type of hazard being modelled, as well as the farm types being considered.

The next step is to combine the damage states into an overall combined impact state (IS) for that event. By allowing for multiple damage states and then combining them into a single overall combined impact state, this framework can be generalised to different hazards and different agricultural land uses. For each overall combined IS, we need to consider how a farm would respond. For most farm types, an event will also impact their ability to operate as usual, and they will have to adjust to a different mode of operation. In many cases, this will result in an increase in expenses, e.g. having to buy extra feed, or an increase in feed or transport prices, and a decrease in income, e.g. due to loss of harvest or inability to get products to market when subject to a transport disruption. These changes are expressed as a vector of scalars, $S_j^{\mathrm{0}}$ (where j is the index of the line item in the farm accounts), that specifies the maximum relative change in each revenue/expenditure line item from the base farm accounts (F_j) due to this behaviour change. We note here that the changes captured in the impact scalars $S_j^{\mathrm{0}}$ can be due to changes in both price and quantity of goods and services. In this way, the model is able to consider the impact of estimated post-event price changes. Developing these line-item-level impact scalars for each combined IS and farm type is a key part of the model for a specific case study region.

At this step, we also need to estimate the impact duration T, which will depend on the overall impact state as well as, potentially, the magnitude of the different hazard impacts. Note here that the impact duration is not the duration of the acute hazard event or even how long the farm's behaviour is affected, but it is the duration of time that the annualised farm accounts are affected. As an example, if there was a disruption that was only 1 week long but it meant that a horticultural farm lost its whole harvest, the impact on the annualised accounts would need to be for a whole year, i.e. T = 1 year.

For our model, we need to calculate the annualised impact on the accounts. For a given year, this means working out the proportion of the year impacted by the event and the mean value of the impact scalar S_j(t) during that time.

The mean value of the impact scalar, during the period of each year that the event impacts, is

$$\begin{array}{}\text{(1)}& {\widehat{S}}_{i,j}={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{t}_i}}\underset{a_i}{\overset{b_i}{\int }}S\left(t\right)\mathrm{d}t,\end{array}$$

where a_i is that start of the event impact in year i ($a_i=min(t^{\mathrm{0}},i-\mathrm{1})$), b_i is the end of the event impact in year i ($b_i=max(t^{\mathrm{0}}+T,i)$), and the duration of the impact in year i is $\mathrm{\Delta }{t}_i=b_i-a_i$. Whether the resulting integral can be evaluated explicitly depends on the form of the recovery function.

The resulting annual farm operating accounts, after adjustments due to the hazard impact, for year i of the simulation are given by

$$\begin{array}{}\text{(2)}& {\widehat{F}}_{i,j}=F_j\times \left[\left({\widehat{S}}_{i,j}\mathrm{\Delta }{t}_i\right)+(\mathrm{1}-\mathrm{\Delta }{t}_i)\right],\end{array}$$

where j is the index for the farm account line items. Note that for all non-impacted periods, S_j=1.

When a farm is impacted by a natural hazard event, in addition to the disruption and changes to business operation, there will often be recovery costs. These can include clean-up costs as well as replacement of damaged assets. We need to determine the total recovery costs C^tot, and for this estimate, we will need farm-specific information, such as building replacement values, effective hectares, and stock numbers. The total recovery costs, C^tot, are then spread out evenly over the impact duration,¹ T, to produce an annualised recovery cost, C_i, for year i:

$$\begin{array}{}\text{(3)}& C_i=\left\{\begin{array}{ll}{C}^{\mathrm{tot}},& \text{if }i-\mathrm{1}<t^{\mathrm{0}}\le i\text{ and }{t}^{\mathrm{0}}+T\le i,\\ \frac{C^{\mathrm{tot}}(i-t^{\mathrm{0}})}{T},& \text{if }i-\mathrm{1}<t^{\mathrm{0}}\le i\text{ and }{t}^{\mathrm{0}}+T\ge i,\\ \frac{C^{\mathrm{tot}}(i-(t^{\mathrm{0}}+T\left)\right)}{T},& \text{if }{t}^{\mathrm{0}}\le i-\mathrm{1}\text{ and }{t}^{\mathrm{0}}+T\le i,\\ \frac{C^{\mathrm{tot}}}{T}& \text{otherwise}.\end{array}\right.\end{array}$$

The output of the model is a time series of the annual revenue, expenses, and costs at the line-item level REC_i,j, for each farm at each location, which is created by appending the recovery costs C_i onto the farm operating accounts ${\widehat{F}}_{i,j}$. The annual net revenue, REC_i, is then ${\sum }_j{\mathrm{REC}}_{i,j}$.

To consider impacts that may be far into the future, the standard economic practice is to discount future cash flows by reducing future cash flows according to a rate of discount over time. This leads to the discounted revenue, expenses, and costs (DREC):

$$\begin{array}{}\text{(4)}& {\mathrm{DREC}}_{i,j}={\displaystyle \frac{{\mathrm{REC}}_{i,j}}{(\mathrm{1}+r{)}^i}},\end{array}$$

where r is the annual discount rate, and i is the year in the future, with i=0 being the point in time for which we are calculating present values. Finally, we calculate the net present value (NPV):

$$\mathrm{NPV}=\sum _{i,j}{\mathrm{DREC}}_{i,j}.$$

2.2 Business recovery

Consistent with a preceding business behaviours model (Brown et al., 2019b), which in that case focused on urban business recovery following an earthquake, businesses are conceived as following a pathway of recovery that is defined by a recovery curve or function across time. However, unlike the Brown et al. (2019b) function, which is applied only to the modification of business outputs or revenues, we require a function that can be equally applied to expenditures. Importantly, depending on the expenditure type, expenditures may potentially be both higher and lower than business-as-usual conditions during the recovery phase. Furthermore, given the range of potential magnitudes for volcanic and other events, it is imperative that the function is flexible to allow for minor as well as severe impacts, as well as short to long recovery durations. Added to this, we prefer a function that is straightforward to evaluate on an annual basis, and it must work in the model for events that happen partway through a year and/or events that have recovery durations completed partway through a year.

Figure 2 Comparison of different recovery functions for an event at t⁰=0.5. In green are the operability functions from the Urban Business Behaviours Model (Urban BBM) (Brown et al., 2019b; Smith et al., 2019) for a range of overall disruption (OD) and delay values. In brown are the impact scalar curves for the new Agricultural Business Behaviours Model (AgriBBM) for a range of values for the impact scalar $S_j^{\mathrm{0}}$, impact duration T, and recovery steepness k.

Given these considerations, we have developed the impact scalar recovery function:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{S}_j\left(t\right)=& S_j^{\mathrm{0}}+\left(\mathrm{1}-S_j^{\mathrm{0}}\right)\exp \left({\displaystyle \frac{k}{T}}[t-(t^{\mathrm{0}}+T\left)\right]\right)\\ & \text{for }{t}^{\mathrm{0}}<t<t^{\mathrm{0}}+T,\end{array}\end{array}$$

where t⁰ is the start of the event in continuous time, e.g. 4.82, T is the impact duration, k is the steepness of the recovery, and $S_j^{\mathrm{0}}$ is the initial maximum value of the relative scalar for farm account line item j. Figure 2 demonstrates the shape of this recovery function for a range of values of T, k, and $S_j^{\mathrm{0}}$, and it includes a comparison with the existing Urban Business Behaviours Model (from Brown et al., 2019b, and including modifications described in Smith et al., 2019).

Whether a recovery curve should follow a concave down-decreasing or concave up-decreasing path has been debated in the literature and will likely depend on the characteristics of the event and sectors impacted (Jonkeren and Giannopoulos, 2014). Given the nature of our case study, with agriculture businesses often facing degraded natural and other capital from which recovery cannot be initiated quickly, as well as with the inability of agriculture businesses to action recovery through relocation, we select a function that captures limited recovery during the initial phases of the impact duration. Nevertheless, the form of the function has some flexibility through recovery steepness parameter, k_e. By setting k < 10, it can represent impacts with more immediate recoveries (almost linear),² and it can represent recoveries that are closer to step functions (shifting back to the base accounts or business as usual only at the end of the impact duration) by setting k≫10. Of course, alternative functional forms could be investigated for other applications.

For this form of the recovery function, we are able to explicitly integrate S_j(t) to calculate the average adjustment to line item j, during the annual accounts period for year i, and thus the adjusted farm accounts. Combining Eqs. (1) and (2), we get

$${\widehat{F}}_{i,j}=F_j\times \left[\underset{a_i}{\overset{b_i}{\int }}{S}_j\left(t\right)\mathrm{d}t+(\mathrm{1}-\mathrm{\Delta }{t}_i)\right].$$

For the chosen recovery function, the integral can be evaluated and is

$$\begin{array}{}\text{(6)}& \underset{a_i}{\overset{b_i}{\int }}{S}_j\left(t\right)\mathrm{d}t={\left[S_j^{\mathrm{0}}t+\left(\mathrm{1}-S_j^{\mathrm{0}}\right){\displaystyle \frac{T}{C}}{e}^{\frac{C}{T}\left(t-\left(t^{\mathrm{0}}+T\right)\right)}\right]}_{t=a_i}^{b_i}.\end{array}$$

2.3 Extending AgriBBM for multiple events

In realistic volcanic scenarios, there will often be more than one event, and the impact of events can be overlapping, e.g. the next event can occur before a farm has recovered from the first event. In our case study, for each hazard scenario, we first split the scenario into specific discrete events e, that occur at times $t_e^{\mathrm{0}}$, and have an associated set of hazard intensity metrics. We now need to extend the single event model to account for these multiple events. In most cases the equations are exactly the same – only an extra subscript, e, added for the event number. However, when events are overlapping, the situation is more complex.

Firstly, for calculating the contributing damage states, DS_e; the overall impact state, IS_e; and the impact duration T_e, we consider each event independently.

For recovery costs, we make the assumption that the cost to repair or replace assets and to clean up will be the same regardless of any preceding events. This means that the total annual recovery costs are $C_i={\sum }_e{C}_{i,e}$, where Eq. (3) becomes

$$\begin{array}{}\text{(7)}& C_{i,e}=\left\{\begin{array}{ll}{C}_e^{\text{tot}},& \text{if }i-\mathrm{1}<t_e^{\mathrm{0}}\le i\text{ and }{t}_e^{\mathrm{0}}+T_e\le i,\\ \frac{C_e^{\text{tot}}\left(i-t_e^{\mathrm{0}}\right)}{T_e},& \text{if }i-\mathrm{1}<t_e^{\mathrm{0}}\le i\text{ and }{t}_e^{\mathrm{0}}+T_e\ge i,\\ \frac{C_e^{\text{tot}}\left(i-\left(t^{\mathrm{0}}+T_e\right)\right)}{T_e},& \text{if }{t}_e^{\mathrm{0}}\le i-\mathrm{1}\text{ and }{t}_e^{\mathrm{0}}+T_e\le i,\\ \frac{C_e^{\text{tot}}}{T_e}& \text{otherwise}.\end{array}\right.\end{array}$$

There is a limitation in our approach for multi-events in that it may overestimate recovery costs, particularly recovery costs that might in reality occur at discrete points in time towards the end of recovery. Take, for example, a farm milking shed damaged during an initial volcanic event. It is possible that a second volcanic event occurs prior to the replacement of this shed being initiated; thus, the replacement cost is only faced for the second event, whereas the model will include replacement costs for the first event. This is a topic for future model refinement, particularly if used for applications where overlapping recovery phases are prevalent in the scenarios considered.

Because the impact scalar vector $S_{j,e}^{\mathrm{0}}$ for each impact state IS_e represents a specific shift in behaviour for the farm, the values for specific line items j cannot be modified independently of each other. This means that combining the impact of overlapping events cannot be done line item by line item; for example, we cannot add the impacts, $S_j\left(t\right)\ne S_{j,e}\left(t\right)+S_{j,e+\mathrm{1}}\left(t\right)$, or take the maximum relative impact of the two event line items. To illustrate this, we can consider the case where there is one event (e=1) causing a farm to lose grazing area and necessitating more bought-in feed to maintain the same stock numbers; here $S_j^{\mathrm{0}}>\mathrm{1}$ for j = feed costs. This can be compared to another event (e=2) where the impact requires the farm to destock drastically, reducing revenue, but also meaning the farm does not need to bring in much feed, if any, so $S_j^{\mathrm{0}}\ll \mathrm{1}$ for j = feed costs.

For our chosen recovery function, the impact scalar at any point in time is the sum of $S_{j,e}\left(t\right)={\sum }_e{S}_{j,e}\left(t\right)$:

$$\begin{array}{rl}{S}_{j,e}\left(t\right)=& S_{j,e}^{\mathrm{0}}+\left(\mathrm{1}-S_{j,e}^{\mathrm{0}}\right)\exp \left[{\displaystyle \frac{k_e}{T_e}}\left(t-\left(t_e^{\mathrm{0}}+T_e\right)\right)\right]\\ & \text{for }{t}_e^{\text{start}}<t<t_e^{\text{end}},\end{array}$$

where $t_e^{\text{start}}$ and $t_e^{\text{end}}$ are set such that only one $S_{j,e}\left(t\right)\ne \mathrm{0}$ at one time, because a farm can only be in one impacted state (described by a single vector of impact state scalars) at a time.

Figure 3 Example plots of how the impact scalar for a specific line item S_j(t) changes through time when two events occur. The dashed black line shows the resulting impact scalar S_j(t). In purple is the impact of event 1, S_j,1(t), with maximum impact scalar $S_{\mathrm{1}}^{\mathrm{0}}$, and an impact duration T₁ indicated by the purple shaded period. In orange is the impact of event 2, S_j,2(t), with maximum impact scalar $S_{\mathrm{2}}^{\mathrm{0}}$, and an impact duration T₂ indicated by the orange shaded period. In panels (a) and (b), event 2 has more impact^∗ than event 1, so if they overlap the farm state changes to the event 2 impact immediately. In panels (c) and (d), event 1 has more impact^∗ than event 2, and if they overlap the farm state stays at the event 1 impact until the end of the impact duration, before switching to the event 2 impact. ^∗ Note that we use recovery time as a proxy for impact magnitude.

When events overlap, we use the impact duration as a proxy for the severity of the event to determine whether the farm stays in the current state or switches to the state of the new event.³ This gives us $t_e^{\text{start}}=max(t_e^{\mathrm{0}}$, $t_{e-y}^{\mathrm{0}}+T_{e-y}\left)\mathit{\right\{}y:T_{e-y}>T_e\mathit{\}}$; that is, if any of the preceding event(s) have a larger impact and the farm is still in recovery from it, this larger event determines the functional state of the farm until its recovery has finished. And $t_e^{\text{end}}=min(t_e^{\mathrm{0}}+T_e$, $t_{e+y}^{\mathrm{0}}\left)\mathit{\right\{}y:T_{e+y}>T_e\mathit{\}}$; that is, when a larger event occurs during recovery, this new larger event supersedes this event (see Fig. 3).

Given the adjustments for when each event e applies ($t_e^{\text{start}}$ to $t_e^{\text{end}}$), we then use that period to work out $a_{i,e}=min(t_e^{\text{start}},i-\mathrm{1})$) and $b_{i,e}=max(t_e^{\text{end}},i)$ for each year i and the mean value of the impact scalar during that period.

$$\begin{array}{}\text{(8)}& {\widehat{S}}_{i,j,e}={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{t}_{i,e}}}\underset{a_i}{\overset{b_i}{\int }}S\left(t\right)\mathrm{d}t,\end{array}$$

where $\mathrm{\Delta }{t}_{i,e}=b_{i,e}-a_{i,e}$.

Finally we calculate the adjusted annual operating accounts, taking into account all events, using

$$\begin{array}{}\text{(9)}& {\widehat{F}}_{i,j}=F_j\times \left[\sum _e\left({\widehat{S}}_{i,j,e}\mathrm{\Delta }{t}_{i,e}\right)+\left(\mathrm{1}-\sum _e\mathrm{\Delta }{t}_{i,e}\right)\right].\end{array}$$

As before, we then create REC_i,j by appending the total recovery costs C_i onto the farm operating accounts ${\widehat{F}}_{i,j}$, and then we can calculate DREC_i,j following Eq. (4).

Figure 4 Overview of the study area: (a) world map showing the location of North Island of Aotearoa / New Zealand; (b) North Island of Aotearoa / New Zealand, highlighting the Taranaki region (StatisticsNZ, 2018a) and Taranaki Mounga; (c) volcanic hazards in the Taranaki region (adapted from Neall and Alloway, 1996) with case study farm locations; (d) key urban areas (StatisticsNZ, 2018b), highways (Land Information New Zealand (LINZ), 2011), and selected rivers (Ministry for the Environment (MFE), 2010) in the Taranaki region; and (e) land use in the Taranaki region (Ministry for the Environment (MFE), 2020). Publisher's remark: please note that the above figure contains disputed territories.

3 Case study

Much of the Taranaki region comprises an extensive ring plain that surrounds Taranaki Mounga (Taranaki volcano) of overlapping volcaniclastic deposit fans (Cronin et al., 2021). A national park, Te Papakura o Taranaki, encircles the volcano out to a radius of approximately 10 km, measured from the volcanic summit. Beyond this distance, the region is dominated by agriculture, with dairy cattle farming being a major land use in the region since around the 1880s (Fig. 4).

The most recent eruptive activity at Taranaki, dated as 1780–1800 CE (Lerner et al., 2019), occurred prior to the development of large-scale agriculture in the region. Volcanism at Taranaki Mounga is characterised by cycles of edifice growth and collapse, with growth periods containing effusive and lava-dome-forming activity interspersed with explosive eruptions (Turner et al., 2011; Lerner et al., 2019; Cronin et al., 2021). Ashfall and the mobilisation of volcanic material in the form of lahars are the two hazards that contribute the most to widespread volcanic risks for agriculture within the region and thus are a focus of this study (see Figs. 4d and 6). Within the volcano's eruptive history, the hazards of lava flows and ballistic projectiles have occurred relatively infrequently (Turner et al., 2011) and have been largely contained within the vicinity of Te Papakura o Taranaki national park, and they are thus not considered in this study. While several large-scale debris avalanches have occurred, caused by the catastrophic collapse of unstable portions of the edifice (Zernack et al., 2012), these events are very infrequent and so are also not considered. We do not explicitly model the risks from pyroclastic density currents (PDCs). The exposure to PDCs is less than that of ashfall and lahar hazards, but PDCs may extend beyond the park boundaries. On the other hand, the impacts of PDCs on farming operations are expected to be similar to lahars, and the footprint of the latter is larger. Thus while initial damage may be due to PDCs, a subsequent lahar would traverse the same area and largely generate the same scale of damage in this case.

Figure 5 Summary of 10 000 volcanic future simulations for Taranaki showing the count of simulations by number of events.

Figure 6 Cumulative ashfall load in the Taranaki region across all 10 000 simulations, on a logarithmic scale.

In terms of cascading failures following a volcanic event, disruptions to water supply and transport networks have been identified as particularly important for dairy cattle farms (Wild et al., 2019). The transport network is dominated by one major highway to each side of the mountain, only three principal access points out of the region, and numerous bridges crossing river channels that have formed in a radial manner around the mountain. Since river channels are likely locations for the flow of lahars, we have focused on nine principal channels in this study (Fig 6). The transport and water supply networks are also subject to risk of damage and disruption from ashfall (Wilson et al., 2014). Farms within the region source water directly from nearby streams, from groundwater, or from community/municipal supply schemes (Wild et al., 2019). The extent to which we expect supplies of water to be disrupted varies according to mode of supply, with those relying on surface water abstraction most at risk of disruption (Porter, 2022).

For this study, we have selected five sites at which dairy cattle farming is located to model the economic consequences of volcanic hazards. The sites were selected to capture a diversity of risk contexts, as well as variation in the land use capability (Lynn et al., 2009). In terms of the latter, we nevertheless note that land use capability or production capacity does not vary significantly among the region's dairy farms, and management choices are likely to be the largest determinant of business-as-usual variation between farm systems.

4 Data and model set-up

4.1 Representative farms

We have not sought to obtain financial data for the real-life dairy farms existing at our selected study locations, as these data would in any case be subject to confidentiality. Instead the approach has been to assign base financial accounts for representative farms at each location. Reference was made to a set of farm financial and environmental data collected for 49 anonymous farms in Taranaki for the 2020–2021 financial year.⁴ With assistance from an agricultural economist involved in the data collection, the accounts for 1 among the 49 farms was assigned to each study location on the basis that the selected accounts described a farm that could feasibly exist at each site. Adjustments were then made to the accounts so that revenue line items reflected long-term average commodity prices rather than the specific prices experienced during the year for which the data were collected.⁵

Table 1 Key attributes and performance indicators for the five case study farms. EBITDA represents earnings before interest, taxes, depreciation, and amortisation. Note that MS represents milk solids.

The five case study farms are representative of dairy farms in the region, ranging from low- to high-input dairy systems, or from system two to five as defined in Hedley et al. (2006). A summary of the key attributes of the case study farms is given in Table 1. There is a wide range in cow performance in terms of milk solids produced per live weight of cows. Farm earnings also vary per hectare, even among farms of the same system type. Analogous results have been found in other studies (see, for example, Moran et al., 2019), highlighting the importance of the unique management and operational choices in determining farm performance.

The five case study locations were allocated water supply types from their locations relative to municipal supply schemes and suitable surface water collection points using the approach from Wild (2016). Properties situated at a distance from municipal supply schemes and suitable surface water collection streams were assigned groundwater collection.

4.2 Volcanic scenarios

A suite of 10 000 simulations (scenarios) of future volcanic activity, each over a 50-year period, were generated for Taranaki. Only eruptions of a significant volume, extending beyond the park boundary, were considered. For each simulation, the first onset time and subsequent repose times (and hence the number of events) were generated from the temporal model described by Damaschke et al. (2018), while the magnitude of each event in terms of volcanic explosivity index (VEI) (Newhall and Self, 1982) is based on the Expert Elicitation model described in Bebbington et al. (2018). The eruption volumes for each event were then obtained by probabilistic interpolation of the VEI, assuming unit differences in VEI to be equivalent to a 10-fold increase in volume.

Of the 10 000 simulations, 57 % contained no volcanic event occurrence, reflecting the underlying probability of a Taranaki volcanic event over the next 50 years (Fig. 5). The majority of the simulations in which volcanic activity occurs contain only one event (30 % of all simulations); however, 10 % contain two events and 2 % three events. The greatest number of volcanic events occurring over a single 50-year simulation is six.

To produce estimates of the ashfall implications of each event, we had available a pre-processed set of 500 potential Taranaki volcanic events, each with an associated spatial ashfall deposition layer for the surrounding landscape (1 km grid). These were produced via runs of the tephra2 (version 2, see https://github.com/geoscience-community-codes/tephra2, last access: 28 October 2024), a tephra dispersion simulation tool with wind velocities sampled from the NOAA Reanalysis database. Each event in our 10 000 simulations was matched to one of these pre-processed ashfall scenarios by sampling from within the ashfall scenarios that had an equivalent volume. To demonstrate that across all simulations ashfall is more likely to be deposited on the eastern side of the mountain due to the prevailing wind directions, Fig. 6 depicts the summation of all ash deposited over all events and simulations for points on or near roads or buildings.

Section 4.4 below provides further information on the processes by which the ashfall information for each event was later converted to metrics of hazard intensity relevant to the damage and fragility functions utilised in the modelling.

4.3 Functions for damage states and impact states

Reflecting the principal risks to farm operations and capital, the combined impact state for dairy cattle farms was calculated by combining the damage states calculated from four separate vulnerability models, each relying on inputs of hazard intensity:⁶ DS^ashfall is the damage state for farms caused by ashfall impacting directly on pastures, crops, farm buildings, and other farm assets; DS^lahar is the damage state for farms caused by a lahar impacting directly on pastures, crops, farm buildings and other assets; DS^transport is the damage state for farms caused by experience of transportation disruption; and DS^water is the damage state for farms caused by experience of loss of water supply services. The derivation of the damage states for each event was as follows:

• DS^ashfall. Here we used fragility functions for pastoral farms from Craig et al. (2021). These fragility functions calculate the probabilities that large pastoral-farming businesses (which include dairy farms) will experience each of the five damage states (refer to Appendix A1), ranging from no disruption (DS0) through to total loss of capabilities (DS4). Thus, given ashfall depth, a damage state could be assigned to each farm and event through probabilistic sampling. Although this model in general allows for ashfall damage state probabilities to be adjusted downwards during times of the year with moderate or low vulnerability, in this case study it was reasoned that the adjustment from Craig et al. (2021) should not be applied given the intensive nature of dairy farming in the region, with highly balanced operations year round.

• DS^lahar. Only two damage states were recognised for a farm with respect to lahars, that is DS0 (no lahar experienced) and DS4 (majority of farm directly impacted by lahar); refer to Appendix A3.

• DS^transport. Fragility functions were constructed following the same structure as the ashfall fragility functions (see Fig. A1b in Appendix A), with the independent variable being the number of days of inaccessibility experienced by the farm (i.e. rather than ashfall thickness). Informed by a previous ex ante study of a major hazard event that had involved identifying key metrics of transport disruption (Smith et al., 2019), the intensity metric of inaccessibility is defined as the inability to travel to any of the three points on a state highway leading into/out of the Taranaki region (see Fig. 4). Altogether, five damage states were recognised (see Appendix A2). While the lower damage state (DS1) describes the situation of transport disruption increasing transportation costs, the higher damage states (DS2, DS3, DS4) also entail loss in the ability to supply output to market and cows dried off and/or culled (DS3, DS4). One aspect to note is that even if the farm location is never completely cut off during an event (inaccessibility duration = 0 d), there is still assigned a high probability for the event of having transport delays and additional costs (DS1), due to ashfall impacts on the wider network. The seasonal nature of dairy farming is incorporated by including a period during which vulnerability is lowered, to reflect when cattle are not being milked.

  Independent of the ashfall-related transport disruptions, we need to consider the transportation disruptions that might be caused by lahars, specifically the loss of bridges making roads impassable. For this case study, we have applied a bespoke adjustment whereby if bridges crossing the principal river channels are required for access to one of the main regional access points and these bridges are impacted by lahar, the probabilities of being in damage states were adjusted upwards. Full details of the calculations are provided in Appendix A2.

• DS^water. Only two damage states were recognised, DS0 (limited disruption) and DS4 (severe disruption). A simple threshold damage function was chosen where DS4 is assigned if the number of days of interrupted water supply for a farm is greater than a specified threshold; otherwise, DS0 is assigned. Details for the calculation are given in Appendix A4.

The combination of damage states experienced by a farm determines the overall combined impact state of that farm. Altogether, 13 unique impact states are recognised. At the extremes, IS0 (no impact) occurs when the damage state is DS0 across all four categories, while IS4 (severe impact, farm ceases operation) occurs when any one of the four categories are DS4. The other 11 impact states represent different combinations of ashfall and transport accessibility damage states. The complete mapping between damage states to impact states is provided in Appendix C.

4.4 Intensity metrics for hazards/cascading failures

4.4.1 Ashfall for DS^ashfall

For each farm and event, the ashfall thickness serves as an input to the fragility function determining DS^ashfall. This was determined by taking the ashfall load from the volcanic scenarios at the closest data point to the centre of the farm. The ashfall load (kPa) was then converted to thickness (mm) using a deposit density of 1000 kg m⁻³, which is commonly used in volcanology studies of eruptions with similar compositions (Barker et al., 2019; Magill et al., 2015; Taddeucci et al., 2011).

4.4.2 Lahar for DS^lahar

There is a very high likelihood of lahars being generated following an eruption event at Taranaki (Procter et al., 2020). The above-described damage functions for lahars require a binary input for each event that describes whether a farm is located directly in the path of a lahar. This information was generated for each event by randomly sampling from probability distributions that relate the likelihood of lahar realisation to (1) the size of an eruption in terms of material deposited and (2) the proximity of the farm to the valleys or catchments in which the ash is principally deposited. For property (1), the data taken for an event from the eruption scenarios are the maximum ash depth as measured around the circumference of Te Papakura o Taranaki national park boundary,⁷ while for property (2) we have measured the radial degrees between the farm location and the point where the maximum ash depth is recorded around the park boundary. The matrices describing the probability distributions are contained in Appendix B1. Note that a unique distribution is provided for each farm, allowing for the probabilities to be set high for those farms located close to the mountain and where the farm, or a large proportion thereof, is located in a valley. Matrices are not provided for Farms 1 and 5, as it was concluded that the lahar risk for these locations is negligible. The assigned probabilities for the three remaining farm locations were based on expert judgement, informed particularly by the risk categories assigned to each location in the lahar hazard mapping (Neall and Alloway, 1996), as well as the results of the modelling of lahars from eruptions of several sizes and weather dispersal mechanisms by Weir et al. (2022).

4.4.3 Inaccessibility for DS^transport

As explained above, inaccessibility is defined as the inability to travel to any of the three points on a state highway leading into/out of the Taranaki region. The road network was split into an existing road segment classification, with these segments ranging in length from 1 m to 28 km, with an average length of 760 m. Segments were assigned the closest ashfall estimate from the spatially explicit eruption simulations. Network analysis was used to check for connectivity between each location and any of the three regional road entry and exit points, assuming no traffic flow on road segments over a threshold of 10 mm of compacted ashfall. Based on recent experience with silt removal in Aotearoa / New Zealand following cyclone Gabrielle (Hawkes Bay Regional Council, 2023), the rate of ash removal is set to 2400 t d⁻¹ across all roads in the region, and that rate was applied evenly over the road network. Connectivity was assessed on a daily time step for each location.

The transport accessibility damage functions also require inputs for each event relating to lahars; in this case it is a binary input for each of the nine principal river valleys (see Fig. 4) specifying whether a lahar occurs, flowing down at least as far as the respective state highway crossings. Provided a lahar reaches this location, it is then assumed the bridge will be damaged. An approach analogous to that used for farm locations was used to assign lahar hazard realisations to the river valleys – refer to Appendix B1 for the assigned probabilities. The probability distribution for Stony River is different from that used for other rivers on the basis that the Stony River valley has a higher likelihood of lahars due to the geography of the valley and as confirmed in geological records (Procter et al., 2020).

4.4.4 Water supply disruption for DS^water

The chosen failure metric for water disruption is the days of lost service. Separate functions were developed for each water supply type to relate ashfall thickness to the estimated duration of water supply disruption, based on the available literature. Interruption for surface and groundwater collection was based on ashfall thickness at the location of the farm. Days of interruption for municipal supply were based on ashfall thickness at the pump power source, water collection point, and the location of the farm. A maximum of 200 d of water disruption was assumed for large events based on the scenario in Porter (2022). Reference can be made to Appendix B2 for further information.

4.5 Impact scalars ($S_j^{\mathrm{0}}$) and impact durations (T)

Impact scalars and impact durations are outlined jointly in this section as there are interrelationships between the parameters selected, with the impact duration being the time over which impact scalars applied. For impact states that are likely of short duration, i.e. less than 1 year, it was considered convenient to generally assign a fixed impact duration, e.g. 1 year, and then apply impact state scalars that would be applicable when considering changes in revenue/expenditure averaged over an entire year subsequent to the event commencement. All remaining impact states were assigned one of five approaches for estimating impact duration – see Table C1 for the complete mapping. In summary, the approaches for estimating impact durations were as follows:

• Ash depth function (DF). Where ashfall is the most impactful hazard realised, the depth of ash at the farm is used to estimate the recovery time by applying a piecewise linear function. The shape of the function, depicted in Appendix C, was informed by data gathered on recovery following the Mt Hudson and Mount St Helens eruptions (Wilson et al., 2011a, b; U.S. Geological Survey, 2024; Dale and Crisafulli, 2018; del Moral and Wood, 1993).

• Inaccessibility period (IP). Where inaccessibility is the dominant impact, the impact duration is either set simply as a defined period (e.g. 0.5 or 1 year) or the length of the inaccessibility experienced by the farm during the milking period.

• Water supply re-establishment (W). Where a farm experiences DS4 for water supply disruption and all other damage states are DS2 or lower, a farm's impact duration is deemed to be controlled by the recovery of water supply. It is considered that by 6 months a farm will have established an alternative water supply option.

• Maximum of ash and water (WDF). Where a farm experiences DS4 for water supply and a high damage state for ashfall, the impact duration is set as the maximum of that calculated from the ash depth function and for water supply re-establishment.

• Maximum of lahar and ash (LDF). For all impact states where a farm has experienced a lahar, the impact duration is set as either that calculated from the ash depth function or for lahar recovery, whichever is larger. An average recovery and rebuild time for lahar is estimated to be 7.5 years, informed by the literature indicating multi-year recovery time frames for pasture and ecosystem recovery from such volcanic disturbances (Crittenden et al., 2003; del Moral and Grishin, 1999; Saputra et al., 2022).

The full set of impact scalars applied in the modelling are available in the Supplement. Note that where a farm experiences IS0, the scalars are all set to 1, meaning no difference from business-as-usual conditions. This can be compared to IS4, where all scalars, except for local government rates payments which are deemed to continue, are set to 0 implying a cease of farm operations.

The derivation of impact scalars for the remaining ashfall-related impact states involved translating the qualitative descriptions of impact states provided by Craig et al. (2021) into best estimates of proportional changes in dairy farm income and expenditure items. In general, as the impact of ash increases, the quantity of feed able to be grown on farms decreases, thus reducing farm revenues and increasing feed costs. At the same time, however, the number of stock remaining on farm also decreases with higher ashfall impacts, providing limited immediate liquidity and reducing the ongoing expenditures necessary for stock maintenance. The same scalars were generally applied across all farms, except in relation to feed and grazing expenditures, as the proportion by which business-as-usual costs for these items change also depends on each farm's starting position regarding the proportion of feed sourced off farm.

For the impact state characterised by only a minor transport disruption (i.e. DS1 for inaccessibility), the scalars for freight and cartage costs were increased. For the impact states made up of more severe transport disruptions (i.e. DS2 and DS3 for inaccessibility), the scalars act to decrease farm revenues, reflecting disruption to production activities.

For impact states that involve a combination of ashfall damage and transport inaccessibility, many of the scalars selected are more severe than the worst case considering ashfall and inaccessibility impacts independently. This reflects that when multiple hazards are experienced, the options available to mitigate impacts reduce. For example, a farmer cannot simply purchase more feed to replace that which cannot be produced on farm due to ashfall damage if transportation to the farm is also disrupted.

4.6 Recovery costs (C_i,e)

In addition to the disruption to normal revenues and expenditures, a farm also potentially faces costs to repair assets and rehabilitate land. For ashfall damage, we assigned these costs to each event according to damage functions that estimate, for each asset type, the proportion of damage experienced based on ashfall depth. Replacement costs for built assets included some margin to cover items such as demolition and debris removal. The main recovery cost for pasture is the cost of ash removal, land rehabilitation, and fencing. Reference can be made to Appendix D for full details. For farms experiencing a lahar, 100 % damage to assets is assumed, with the same replacement and land rehabilitation costs applied as for ashfall.

A further and potentially significant cost faced by farms following an event is the cost to restock cattle. Implicit to the setting of the impact scalars that modify farm revenue and expenditure line items (see Sect. 4.5) was the assigned proportion of stock that can remain on each farm under each impact state (see Table C1 in Appendix C). When a farm experiences a damage state of DS2 or below for transport, it was assumed that half of the value of the stock that must be removed from the farm can be recaptured by stock sales. However, for the higher transport damage states, it was assumed that farms have no opportunity to recoup any value of the cattle that cannot be supported on-farm given the transport difficulties in getting stock to market.

5 Procedure

We ran AgriBBM for each of the 10 000 volcanic scenarios and for the modelled dairy farms (i.e. Farms 1–5). We also tested locating each of the five modelled farms at each of the five study locations. Note that it is not necessarily the case that each of the farms chosen as representative farms could be located at any of the study locations, given that site characteristics were taken into consideration when selecting the representative farm for each location. Nevertheless, it was informative to test running the model for each farm at each location because it allowed for separating the implications of farm management (e.g. proportion of imported feed) from the implications purely related to farm location (e.g. proximity to vent). The mode of water supply available to a farm (e.g. surface water versus groundwater collection) was held fixed at a location. We also fixed the random numbers generated for fragility function sampling by event and location in order to best enable comparisons across farms placed at the same location.

Figure 7 Annual net revenue (NZD thousands) for five farms at five locations for a selected scenario which included three eruption events. Icons represent the incidence and magnitude of volcanic impacts.

To illustrate the results from each model run, we show 50 years of annual net revenues for an example scenario in Fig. 7. This scenario had three eruptions within the 50-year simulation, occurring at 20, 22, and 42 years. The depth of ashfall and the realisation of lahars varied across the five locations. The figure depicts separately the adjusted annual farm operating profit and the total annual recovery costs, with the summation being the net revenue each year.

6 Results and discussion

In this section, we summarise the results across the 10 000 ashfall simulations and discuss some of the findings and priorities for future work. Although it is possible to draw out detailed (line-level) results for each farm, we concentrate on the damage states and the overall revenues, expenditures, and costs experienced by the farms under the various simulations.

Figure 8 The proportion of damage states among single event simulations for Locations A–E: (a) ashfall damage states, (b) lahar damage states, (c) transport damage states, and (d) water supply damage states.

6.1 Damage states

The impact state of farms was determined by combining the four damage states calculated for ashfall, lahar, transport, and water supply. Figure 8 describes the distribution of damage states at each of the five farm locations, across the events. For simplicity, in undertaking this and the subsequent analysis in Sect. 6.2, we have only considered scenarios where there is a single event during the 50-year simulation.⁸ Note that the assignment of the damage state was by location, and was independent of the farm.

The results for Locations B and D are clearly dominated by the frequent incidence of lahar, with lahars associated with 37 % and 66 % of all events at these respective locations. Given that these sites are the closest to the vent and located on the side of the mountain where ash deposition is more prevalent (especially at Location D), it is not surprising that these locations also experienced the most severe damage states in terms of ashfall and transportation.

Although Location E is on higher ground with negligible lahar risk, being located to the east of the mountain, it does experience quite significant ashfall and transportation risks, with 35 % of events producing a damage sate of DS2 or higher for ashfall and 50 % DS2 or higher for transport. Located at some 37 km from the mountain and also to the northwest of the mountain, Location C experiences overall the least severe ashfall damage states. Despite non-zero probabilities for lahars, none of the simulations produced lahars affecting Location C due to the lighter ashfall and ashfall-driven hazards on the western side of the mountain. Nevertheless, with many kilometres of road and channels to cross in either direction to gain access to or from the site, there were significant risks for transportation disruption. Additionally, for Location A, the loss of water supply appeared to be a significant risk, with 42 % of events returning a damage state of DS4. The provision of more resilient back-up water supply options (e.g. groundwater) appears highly worthwhile at locations similar to Location A, as it would significantly improve the risk profile for these locations.

With lahars shown to be such an important component of the risk profile for farms located in close proximity to the vent and in river valleys, it is an important topic for further hazard research in Taranaki. The way in which lahar realisations were allocated to our scenarios has been based on rather broad assumptions and probability matrices, albeit based on some prior research of the hazard. With more developed hazard information, we could expect some more nuanced findings for the different farms and locations, although this might also require some modifications to our modelling approach. One aspect in which the approach could be more nuanced is in relation to the timing of lahars, with these potentially being treated as separate events. It might also be possible to define different recovery times depending on the material deposited from (and potentially the velocity of) a lahar. Another aspect worth considering is the community or governance decisions that might be put in place in response to lahar risk, such as evacuations, and the resulting economic implications for farms.

Figure 9 Total economic impact of single (non-overlapping) volcanic events as measured by total cumulative change in net annual revenue (thousands of NZD  per hectare), across five farms (1–5), at five locations (A–E).

6.2 Location versus farm

To help with understanding the importance of location relative to farm characteristics in determining economic risk, the box-and-whisker plots in Fig. 9 describe the distributions of economic impacts across events, with the model run for all farms at all locations. Here the economic impact of an event is derived as the (non-discounted) difference between the sum of the annual net revenues (∑_iREC_i), as calculated for the scenario in which the subject event occurs, and the scenario where there is no volcanic event. To help adjust for the influence of farm size on the resulting economic impacts, the economic impacts are reported per hectare.

It is apparent from Fig. 9 that variation among farms has a small influence on economic impacts compared to farm location. For Locations A, C, and E, there is the classic long-tail distribution of event impacts. This distribution shape with a high proportion of zero or close-to-zero losses occurred despite all scenarios having volcanic events. In particular, for Locations A and C, all farms exhibit a median event impact of zero loss. At location E, the median event impact was slightly higher, reflecting the location's greater ashfall and transportation risks, but still less than NZD 2000 per hectare of loss across all farms. Each of these three sites nevertheless has the potential to experience substantially higher economic impacts with the rare but significant events represented by the outliers in Fig. 9 of over NZD 10 000 per hectare.

At Locations B and D, the incidence of lahars causes distributions with a significantly different profile. At Location D, with lahars occurring across more than half of the events and with the calculated economic impacts of such events all being the same, we see both the median and 75th percentile economic impacts falling at the ends of the respective distributions for each farm. Location B also has a large proportion of scenarios with severe impacts, and indeed the 75th percentile economic impacts are the same as calculated for Location D (i.e. ranging from NZD 28 000 per hectare to NZD 38 000 per hectare across the five farms).

While there is generally little variation among the farms when placed at the same location, we can note that at Location A the 75th percentile economic impact for Farms 2, 4, and 5 is less severe than that of the other farms. Given that we have already controlled for farm size, these lower impacts partly reflect the lower stocking rates on these farms, meaning there are lower costs for asset and stock replacements. Farm 4 is also shown to experience lesser impacts at other locations, with, for example, the 75th percentile economic impact for Farm 4 the least severe at Locations B and D. To a large extent, this reflects that Farm 4, being the largest farm, has the most moderate revenue per hectare (i.e. the lowest per hectare EBITDA), and so there is simply less to lose per hectare when the farm is out of operation, as well as low costs per hectare to recover assets and replace stock.

In this analysis, we have accounted for farm heterogeneity in terms of base characteristics of the farm before an event, e.g. stock numbers, proportion of feed imported, cost and revenue structures. In terms of recovery pathways following an event, we have simply incorporated the best approximation of an average recovery curve for farms and assumed a farm always returns to its prior state. However, given the observation that management decisions have a significant influence on farm performance under business-as-usual conditions, it is reasonable to propose that management decisions might also have an influence on the economic impact of an event. Thus, further research could be devoted to identifying the range of potential pre-emptive actions and responses to an event and quantifying the economic implications of these different management decisions. Although there might be little ability for farmers to avoid damages to a farm's natural and physical capital, there potentially are choices around the management of stock (e.g. percentage of stored feed, feed shelters, and, if possible, pre-emptive evacuation when volcanic unrest is high), investment in land rehabilitation, and selection of alternative types of management practices following an event to accelerate recovery (including even selection of alternative farm types).

Figure 10 Net present value of the economic impacts over the next 50 years for five case study farms, as simulated under 10 000 scenarios and with a discount rate of 6 % per annum. Note that in 57 % of scenarios, there are no volcanic events; thus, the impact is NZD 0.

6.2.1 Economic impacts of volcanic hazards over the next 50 years

Figure 10 presents the distribution of economic impacts across the scenarios for the entire 50-year simulation period. Here the economic impact of volcanic hazards for a scenario is calculated as the difference between the net present value (NPV) of that scenario and the NPV of a scenario with no volcanic events. Unlike the previous Fig. 9, where all five farm types were shown at each of the five locations, Fig. 10 displays each farm only in its intended location. Note further that all 10 000 simulations are included in this figure, encompassing multi-event scenarios as well as the 57 % of scenarios that do not involve any eruptions. Given the latter observation, the 25th percentile and 50th percentile of losses fall to zero for all farms. A further preliminary observation is that discounting does itself help to spread the impacts under the scenarios, with higher losses calculated for scenarios having events that occur in the near future compared to scenarios where events occur in the distant future.

Notably in Fig. 10, farms located in the two lahar-affected areas show a longer tail of severe impacts across the simulated volcanic events. These two farms also have higher mean NPV losses, with Farm 4 at Location D a mean loss of NZD 6080 per hectare and Farm 2 at location B a mean loss of NZD 3570 per hectare. In contrast, Farm 3 at Location C, which is not in a lahar-affected area, shows a lower mean NPV loss of NZD 540 per hectare. Despite not being affected by lahars, Farm 5 at Location E experiences a mean impact due to its high ashfall exposure of NZD 1370 per hectare.

When it is considered that the average value of dairy cattle farms in Taranaki is NZD 44 000 per hectare (Colliers International, 2023; Real Estate Institute of New Zealand, 2024), the risk of economic losses to dairy farms in Taranaki is shown to be far from negligible. The findings underscore the varying degrees of exposure across different farm locations and the significant role of volcanic event proximity and local terrain in determining the extent of economic risks.

Figure 11 Probability density functions of the economic losses from volcanic events over the next 50 years for five case study farms. Note that (1) VaR represents the value at risk for the 95th confidence level and (2) NPV loss represents the net present value of total losses in net annual revenues when applying a discount rate of 6 % per annum.

6.3 Value at risk

The value at risk (VaR) statistic is a useful risk management tool to assist in conceptualising risk. Essentially, it defines the worst potential losses that can be expected, over a defined time period, for a given confidence level. The ridgeline plots of probability densities in Fig. 11 illustrate the placement of the VaR for a 95 % confidence level. Note that to enhance the visibility of the larger NPV losses, the y axis of Fig. 11 has been transformed using a fourth-root transformation. This magnifies the lower end of the distributions, making less common high-impact scenario densities more visible.

Using the 95 % confidence level VaR, Farm 2 and Farm 4 exhibit the highest value at risk from volcanic events, primarily due to lahar impacts. Among Farms 1, 3, and 5, Farm 3 has the lowest value at risk. The value at risk across the five farms ranges from NZD 5200 per hectare to NZD 35 700 per hectare, with non-lahar-affected farms reaching NZD 8100 per hectare and lahar-affected farms ranging from NZD 21 300 per hectare up to NZD 35 700 per hectare.

The ridgeline plot (Fig. 11) highlights that the longest distribution tails correspond to the most impacted locations. Notably, the simulation with the maximum financial impact reached NZD 93 000 per hectare for the most heavily impacted Farm 4. For non-lahar-affected farms, the maximum was NZD 29 700 per hectare at Farm 5. The extreme cases will correspond to scenarios with not only high impact and near-term events, but also incidences of multiple volcanic events across the 50-year horizon.

Figure 12 Value at risk (VaR) curves for five case study farms derived from simulated NPV losses over the next 50 years.

The VaR for alternative confidence levels can also be calculated, which is informative for decision-makers and communities who will have a range of risk preferences. In Fig. 12, we plot the VaR for increasing confidence levels. Due to the predominance of no-eruption scenarios, the VaR for all five farms remains zero until nearly 60 % confidence. After this point, it begins to rise, particularly for Farm 4 but also for Farm 2. As in Fig. 12, Farm 4, the most impacted farm, has a 90 % confidence level VaR of NZD 22 100 per hectare, increasing to NZD 71 600 per hectare at 99.9 % confidence. For the non-lahar-affected Farm 5, the 90 % confidence VaR is NZD 4800 per hectare, increasing to NZD 22 200 per hectare at 99.9 % confidence. The results for Farm 3 indicate that farms located away from the main ashfall areas to the east of the volcano and with low lahar risk experience losses only in the worst 10 %–15 % of scenarios (VaR confidence over 85 %).

7 Conclusion

The quantification of economic risks to farms or other businesses from volcanic hazards is a multi-disciplinary, system-oriented exercise. We are required to pass appropriate information and research findings from the underlying physical volcanic processes through models of hazard generation and models capturing cascading failures of support systems, prior to being able to even consider the dynamics of economic impacts and recovery. This study has provided a successful demonstration of such an end-to-end endeavour. It has also produced an agricultural business model that is sufficiently general to be applied in other similar contexts.

We have shown that dairy cattle farms in the Taranaki region are generally sufficiently similar such that variations in farm management will likely have a small influence on economic risk when compared to the variation in risk attributed simply to spatial variation in hazard exposure. For the three case study farms that had negligible lahar exposure, the modelling results indicate that, with a level of 90 % confidence, the economic impacts from volcanic events over the next 50 years will not exceed approximately 10 % of the property value of those farms. By comparison, in the case of the farm with the most significant lahar and ashfall exposure, the results indicate that we cannot be more than around 80 % confident that the impacts over the next 50 years will not exceed around one-quarter of the farm's value, and we also cannot be more than around 90 % confident that the impacts will not exceed approximately half of the farm's value. These results indicate that, provided the region's dairy sector has sufficient access to risk information, we should anticipate volcanic risk as having a significant influence on investment choices and other management decisions that shape the future spatial evolution of the sector.

The modelling pipeline and risk metrics demonstrated in this paper have been developed to be sufficiently generalisable to be able to be used to investigate a number of important research questions. A priority for future research will be to use these tools to assess the impact of potential mitigation and adaptation strategies intended to reduce the risk from volcanic hazards, as well as to improve the preparedness and resilience of farms. For example, one potential pre-emptive measure that could be investigated is land use change for areas of high exposure (e.g. to arable, cropping) and whether this produces any benefits in terms of reduced economic and animal welfare risks. This model could also be used to consider the appropriateness of ex ante self-insurance strategies such as diversification of income sources, increasing the share of capital dedicated to off-farm activities, or geographically spreading the risk by investing in farms outside of the risky zones (Choumert-Nkolo et al., 2021). In these regards the value-at-risk metric used in this study for a self-insurance scenario(s) would be estimated and compared against the no mitigation scenario as part of an evaluation of the merits of mitigation strategies.

Agricultural activities are highly seasonal in nature, and the impact of an event on a farm (and on the farm accounts) can be highly dependent on the time of year. An important limitation of this work is the reliance on annualised farm accounts. The development of seasonal line-item level farm accounts would enable a more realistic consideration of the seasonal impact of events. This would be important in particular for modelling hazard scenarios that are themselves seasonal in nature, e.g. severe weather events, or for farm types such as horticulture, which are more seasonal. Despite the limitations, in this case study, we have considered some of the primary seasonal effects for dairy farms, specifically we consider the overlap of any event with the milking season to determine both the transport damage state (and final impact state) and to calculate the duration of the impact on the farm accounts. In future applications of this modelling pipeline, careful consideration needs to be given to the way in which impacts will vary depending on the time of year the event occurs.

There are several opportunities to refine the model and analysis in future, e.g. through integration of more developed lahar hazard information, the ability to consider compounding events, and more nuanced consideration of impact duration and recovery costs for overlapping events. Additionally, the current model focuses on impacts to farms following the realisation of discrete hazard events. Taking a wider view of what may constitute an event and also recognising that life-safety will be a high priority for many government and non-government agents, it would be useful to also consider the disruption to farms caused by evacuation orders and decisions. Evacuations may occur when there are signs of increased volcanic activity or when lahar risk is heightened, and evacuations are not necessarily followed by hazard realisation. Overall, recognising there are areas for future extension and refinement; the findings presented in this paper are nevertheless among the most informative that are available to assist farmers and other decision-makers in understanding the risks to dairy cattle farming from volcanic hazards in the Taranaki region.

Appendix A: Fragility functions

A1 Ashfall

Ashfall damage states are described as the following:

DS0	No disruption
DS1	Some disruption
DS2	Moderate disruption
DS3	High disruption
DS4	Total loss of capabilities

We use the piecewise linear fragility function parameters for large pastoral farms from Craig et al. (2021), with the additional assumption that for ashfall depths greater than 1000 mm we assume the damage state is always DS4. These are shown in Fig. A1a, and the values for the points on the graph are given in Table A1.

Figure A1 Probability that damage state is achieved or exceeded due to the following: (a) ashfall thickness (DS^ashfall) and (b) days of road inaccessibility (DS^transport). Note here that even for a situation with no time inaccessible (inaccessibility = 0 d), there is still a high probability of being in DS1, due to wider network disruptions.

Table A1 Boundary values for the piecewise smooth ashfall fragility functions.

A2 Transport

The five transport damage states are described as the following:

DS0	no impact;
DS1	some disruption – additional costs faced by farms;
DS2	moderate disruption – some output from farms cannot get to market;
DS3	high disruption – feed supplies depleted, some stock culled, loss of output to market, cows dried off;
DS4	severe disruption – farms cease operations, cows culled.

We have developed piecewise linear fragility functions for dairy farms that depend on the number of days that the farm is completely inaccessible by road due to ashfall impacts on the roads. These fragility functions are shown in Fig. A1b, and the values in the plots are given in Table A2. It is important to note here that even for a situation with no time inaccessible (inaccessibility = 0 d), we assume that there is still a high probability of being in DS1, due to wider network disruptions. This is different to many other fragility functions in the literature.

The time of year of the event will have a significant effect on how impacted a farm will be. For dairy farms, the key consideration is whether the transport disruption falls during (or extends into) the milking season. To account for this, the probability for each damage state calculated from Fig. A1b is adjusted downwards for transport disruptions that extend over the lower vulnerability period (see Table A2 for details).

In addition to the disruption caused by ashfall on the roads, in this region, we also need to take into consideration the transport disruptions caused by bridges on key roads being impacted by lahars. Here we consider bridges on the state highway network that cross nine principal river channels. If any one or more of the bridges en route to any regional access point are damaged, we increase the probability of being in DS1. Then, if at least one bridge on each side of the farm is damaged, this prevents access out of the region, and we increase the likelihood of being in DS2 and of being in DS3. The probabilities vary by farm location and are given in Table A3.

The final transport damage state probabilities are the maximum of the inaccessibility due to ashfall fragility function calculation (Table A1) and the lahar bridge impact (Table A3).

A3 Lahar

The two damage states for lahar impact are the following:

DS0	limited or no impact from lahar,
DS4	majority of farm directly impacted from lahar.

Here we simply set the damage state to DS0 if the farm is not directly hit by a lahar and to DS4 if it is.

A4 Water supply

We only consider two damage states for loss of water supply:

DS0	limited or no disruption,
DS4	severe disruption – no water available for long enough that the farm has to dry off or sell cows and cease operation.

Here we applied a simple threshold function, where if there is no water supply for 30 d or more the farm is in DS4 or otherwise DS0.

Table A2 Boundary values for the piecewise smooth ashfall fragility functions, as well as timing and adjustment parameters for low vulnerability period(s).

Notes: ^∗ This is the maximum adjustment. The adjustment becomes less when only a portion of the inaccessibility phase occurs in the lower vulnerability period.

Table A3 Transport damage state probabilities due to lahar impacts on bridges on the state highway network, for different locations.

Notes: ^∗ Only bridges on state highways crossing nine principal river channels are considered.

Appendix B: Hazard realisation

B1 Lahar realisation

Given the distribution of ashfall from an event, we used the probabilities specified in Tables B1 and B2 to, respectively, assign the realisation of lahar to each farm and principal river valley.

B2 Water supply failure

The length of time for water supply failure was determined from ashfall thickness at the water supply nodes relevant to the water supply type: municipal, surface water collection, or groundwater collection as shown in Fig. B1. We assigned supply types to farms based on location relative to municipal supply schemes and suitable surface water sources, following Wild (2016). Damage functions were derived from information on municipal supply impacts in Wild (2016) and from farm impacts described in Porter (2022), Wilson et al. (2014), Thompson et al. (2017), and Neild et al. (1998). Specifically, the function for estimating days without water supply for municipal systems depended on the ashfall thickness at the farm, municipal supply intake, and local power transformer station. Meanwhile, the functions for estimating days without surface and groundwater supply were based on the ashfall thickness at the farm.

Table B1 Probability that farm will be impacted by lahar.

Table B2 Probability that river will be impacted by lahar at a point of state highway crossing.

Figure B1 Water availability outage length in days as a function of ashfall thickness for (a) municipal water supply as the maximum delay across delays at the water collection point, farm, and local power source (power station) and (b) ground and surface water sources dependent on ashfall thickness at the farm.

Appendix C: Impact states and impact duration

The mapping between damage states, impact states, and methods used to determine impact duration is provided in Table C1. For each impact state, the share of stock remaining on site is also provided in Table C1. Where stock were removed, it was assumed only 50 % of their sale value could be realised because they are selling under duress. Note that these shares were important in deriving the impact scalars for each impact state. Due to the size of the matrices, the full set of impact scalars is provided in the Supplement. The function that was used to approximate impact duration given ashfall depth (DF) is described in Fig. C1. Impact duration was based on descriptions and data on recovery from previous eruptions (Wilson et al., 2011a, b; U.S. Geological Survey, 2024; Dale and Crisafulli, 2018; del Moral and Wood, 1993). In the case of experiencing a lahar, it is possible that animal mortality will not initially be 100 %. It is nevertheless assumed that no stock sales can be realised due to the likely concurrent impacts on the transport network.

Figure C1 Event impact duration as a function of ashfall thickness for forest and non-forest land uses.

Table C1 Mapping between damage states, impact states, impact duration, and stock outcomes.

Notes: 0 yr = 0 years, 0.5 yr = 0.5 years, 1 yr = 1 year, DF = ash depth function, IP = inaccessibility period, W = water supply re-establishment, WDF = maximum of ash and water, LDF = maximum of lahar and ash, ifp = proportion of farm feed that is normally imported, n/a = not applicable.

Appendix D: Recovery/clean-up costs

Asset replacement costs

The damage functions for farm assets are described in Fig. D1, while the parameters establishing the replacement or rehabilitation costs for assets are provided in Table D1. Replacement values are intended to include some associated costs. The insurance replacement calculator used for the farm dwelling includes the cost of structures, special features, labour, professional fees (e.g. architects, council consents), demolition costs, and removal of any debris (Vero Insurance New Zealand, 2024). Any further associated costs have not been included due to lack of information.

(Vero Insurance New Zealand, 2024)(Ministry for Primary Industries, 2016)

Table D1 Recovery costs for farm assets.

Figure D1 Damage ratios used for calculating losses as a proportion of asset value. Damage ratios were a function of ashfall thickness (mm) at the farm determining losses for buildings, fences, machinery, and sheds.