Asef Nazari
email: asef.nazari@csiro.au

Andreas Ernst
email: andreas.ernst@csiro.au

Simon Dunstall
email: simon.dunstall@csiro.au

Brett Bryan
email: brett.bryan@csiro.au

Jeff Connor
email: jeff.connor@csiro.au

Martin Nolan
email: martin.nolan@csiro.au

Florian Stock
email: florian.stock@gmex.de

Combined Aggregation and Column Generation for Land-Use Trade-Off Optimisation

Keywords: aggregation-disaggregation, column generation, clustering, land use allocation adfa, p. 1, 2011

In this paper we developed a combination of aggregationdisaggregation technique with the concept of column generation to solve a large scale LP problem originating from land use management in the Australian agricultural sector. The problem is to optimally allocate the most profitable land use activities including agriculture, carbon sequestration, environmental planting, bio-fuel, bio-energy, etc., and is constrained to satisfy some food demand considerations and expansion policies for each year from 2013 to 2050. In this research we produce a higher resolution solution by dividing Australia's agricultural areas into square kilometer cells, which leads to more than thirteen million cells to be assigned, totally or partially, to different activities. By accepting a scenario on agricultural products' return, carbon related activities, future energy prices, water availability, global climate change, etc. a linear programming problem is composed for each year. However, even by using a state of the art commercial LP solver it takes a long time to find an optimal solution for one year. Therefore, it is almost impossible to think about simultaneous scenarios to be incorporated, as the corresponding model will become even larger. Based on the properties of the problem, such as similar economical and geographical properties of nearby land parcels, the combination of clustering ideas with column generation to decompose the large problem into smaller sub-problems yields a computationally efficient algorithm for the large scale problem.

Introduction

Predicting land use change in Australian agriculture in the context of increasing energy prices stimulating bio-fuels and bio-energy land uses, and a carbon policy with possibilities for increasing carbon price over time requires solving a large scale linear programming problem. The model covers the domain of southern and eastern Austral-ian land currently in intensive agricultural use on a one square kilometer grid cell resolution over 813,000 square kilometers resulting in more than thirteen million square cells, including active and inactive cells, as shown in Fig. 1. Starting with the agricultural land present in 2013, the model finds the optimal land use each year forward in time until 2050 where optimality is defined as maximizing profit and social welfare (the sum of profit and consumers' surplus). In the highest resolution with land parcels of size one square kilometer, for each year the corresponding LP has 7,313,847 continuous variables and 814,811 constraints. The focus is on the change from current agricultural production to alternative land uses such as carbon plantings, environmental plantings, bio-energy or bio-fuels. Food prices are computed endogenously following a maximum welfare approach whereas costs and revenues for nonagricultural commodities depend on scenarios and exogenous modeling. The platform is built to model a range of scenarios involving alternate assumptions about global climate change, world carbon emissions trajectories, emissions limits and prices of carbon credits, price trajectories for energy, world food demand, supply and price trajectories, and agricultural productivity growth [START_REF] Bryan | Supply of carbon sequestration and biodiversity services from Australia's agricultural land under global change[END_REF].

Fig. 1. Area for land use trade optimisation

Efficiently solving a large scale LP problem is the main difficulty that we address in this paper. At the present setting, finding an optimal solution for each year takes a long time even using commercial software packages. Also, considering that the model is build based on one set of fixed scenarios, the further expansion of the project to consider simultaneous scenarios for the sake of stochastic optimisation is almost impossible at this stage. Therefore, our first aim is to develop an algorithm which can solve the large LP pertaining to a particular year in a shorter amount of time.

One important fact in terms of agricultural land allocation for different products is that nearby land parcels have similar economical and geographical properties. They indeed have similar proximity to water resources, and main roads leading to similar transportation and production costs. In reality it is easy to detect that huge parcels of land are allocated to the same agricultural commodity while passing large farms in countryside. By drawing inspiration from this fact, instead of considering cells of size one square kilometer, the whole Australia is divided into big chunks of land (lower resolution) called clusters. The economical properties of each cluster are calculated based upon cells in each cluster. This type of aggregation technique in optimisation is a handy tool to create a set of smaller problems out of a large problem. However, the smaller problem size comes at the cost of lower accuracy. We remedy this by using an iterative disaggregation approach. The smaller problem is gradually made a better approximation of the original problem by breaking down big clusters into smaller clusters, and add new clusters as new column to the small problem.

There is an increasing level of attention to use optimisation techniques in agricultural and other land use management planning context. As an example, in [START_REF] Weintraub | A Hierarchical Approach to Forest Planning[END_REF] a hierarchical approach is presented for large-scale forest planning. The algorithm is based on solving an aggregate problem, which is of moderate size. Another example, in [START_REF] Seppelt | Identifying trade-offs between ecosystem services, land use, and biodiversity: a plea for combining scenario analysis and optimization on different spatial scales[END_REF] authors argue the usage of optimisation techniques in combination with scenario analysis can provide efficient land use management options for sustainable land use from global to sub-global scales. In terms of water resource management, in [START_REF] Gaddis | Spatial Optimization of Best Management Practices to Attain Water Quality Targets[END_REF] a spatial optimisation techniques implemented among four diffuse source pathways in a mixed-use watershed to maximize total reduction of phosphorus loading to streams while minimizing associated costs. An interesting utilization of a multi-objective optimisation technique is reported for identifying optimum land management adaptations to climate change [START_REF] Klein | Adapting agricultural land management to climate change: a regional multi-objective optimization approach[END_REF].

The rest of the paper is arranged as follows. In the next sections, we describe the model in more details. Then a short description of column generation and aggregation techniques is provided. The new algorithm is introduced at the next section, and numerical results and the conclusion will conclude the paper.

The Model

In order to insert demand-production equilibrium, a segmentation of the domain is necessary to approximate the nonlinear relationship between the demand and production using piecewise linear functions. Also, an appropriate model needs to consider satisfaction of the Australian agricultural food demand by imposing supplyproduction and supply-demand constraints. As we should allocate every cell, partially or totally, to the activities, a set of land-use constraints are introduced. Furthermore, in case there are some particular expansion plans for particular activities, some new expansion constraints are introduced to fulfill existing capacity for each activity. The complete model of the land-use trade-off optimization (LUTO) project is represented by equation [START_REF] Dantzig | Decomposition principle for linear programs[END_REF]. Some necessary and technical information of the model is represented in Table 1.

With resolution set to one, which is the highest resolution with land parcels of size one km 2 , there are |R|=812,383 regional active cells. Also, there exist |S|=100 price and demand segmentation, |J|=9 activities and |F|=24 food commodities.

Constants

𝑄𝑄 𝑠𝑠𝑠𝑠 Quantity of activity j that can be produced at cell r.

𝛿𝛿 𝑟𝑟𝑠𝑠

Annualized economic return of activity j at r.

𝐷𝐷 𝑠𝑠𝑠𝑠

Demand for commodity j at segment s, ∑ 𝐷𝐷 𝑠𝑠𝑠𝑠 = 𝐷𝐷 𝑠𝑠 𝑠𝑠 .

𝐶𝐶 𝑟𝑟𝑠𝑠

Production cost of commodity j at region r.

𝑃𝑃 𝑠𝑠𝑠𝑠

Price of commodity j in segment s at demand level 𝑦𝑦 𝑠𝑠𝑠𝑠 .

Table 1. Components of the model

Column Generation

Column generation is a widely used technique to solve large scale linear and integer programming problems starting from pioneering publications [START_REF] Dantzig | Decomposition principle for linear programs[END_REF] and [START_REF] Gilmore | A linear programming approach to the cutting-stock problem[END_REF]. The technique is frequently utilized when the number of variables is much larger that the number of constraints. In such large scale problems, the vast majority of the variables are zero at optimality, hence the fundamental concept underlying column generation is to solve a smaller problem instead of the original LP by considering a subset of columns (variables). The smaller problem is referred as restricted master problem (RMP), and new columns are added as required. The generation of new columns is accomplished by solving another problem called pricing sub-problem following each optimisation of the RMP. In the pricing stage of each iteration, the column(s) with minimum (maximum, depending on the objective) reduced cost(s) are added to the RMP. The optimality is achieved if it is impossible to add a new column to the RMP. For comprehensive surveys on column generation interested readers can consult [START_REF] Desaulniers | Column Generation[END_REF] and [START_REF] Lübbecke | Selected topics in column generation[END_REF].

For a formal explanation of the column generation algorithm consider the linear programming [START_REF] Dantzig | Decomposition principle for linear programs[END_REF] where the 𝒑𝒑 𝑠𝑠 is the 𝑗𝑗 𝑡𝑡ℎ column of the coefficient matrix, and the number of variables is extremely greater than the number of constraints. Assume that an initial basic feasible solution, 𝒙𝒙 𝑒𝑒 , is available, with associated basis matrix 𝐸𝐸, and cost coefficient 𝒄𝒄 𝑒𝑒 . The simplex multipliers associated with this basis could be calculated as 𝜋𝜋 = 𝒄𝒄 𝑒𝑒 𝐸𝐸 -1 and are always available at each iteration of the simplex algorithm. To improve the basic feasible solution we "price out" all columns corresponding to non-basic variables by forming their corresponding reduced cost 𝑃𝑃 𝚥𝚥 � = 𝑃𝑃 𝑠𝑠 -𝜋𝜋𝒑𝒑 𝑠𝑠 . If max 𝑃𝑃 𝚥𝚥 � = 𝑃𝑃 𝑠𝑠 � > 0, then considering non-degeneracy, the current solution may be improved by introducing 𝑥𝑥 𝑠𝑠 into the basis via a pivot transformation.

Aggregation Techniques

An important issue in obtaining an optimal solution of large scale optimisation problems is the trade-off between the level of details and the ease of solving the model in an acceptable amount of time. Aggregation-disaggregation techniques provide some methodologies for handling large optimisation problems by combining data, or using aggregated problems which are reduced in size. One main approach to construct an aggregated problem for a large scale LP is by partitioning the variables and forming corresponding columns by weighted average of columns in each partition. A feasible solution of the original problem is then achieved by applying a special transformation, called disaggregation, to an optimal solution of the aggregated problem. The aggregation error calculated as the difference between the original optimal objective value and the optimal value for the disaggregated solution guides the algorithm towards the optimality in an iterative scheme. In this section, a concise description of aggregation inspired from [START_REF] Litvinchev | Aggregation in Large-Scale Optimization[END_REF] is provided. Consider the original linear programming problem in the form of

(𝐿𝐿𝑃𝑃) � 𝑧𝑧 * = max 𝒄𝒄𝒙𝒙 𝑠𝑠. 𝑡𝑡. 𝐴𝐴𝒙𝒙 ≤ 𝒃𝒃 𝒙𝒙 ≥ 𝟎𝟎 (2)
A column aggregation is explained here, which means only variables in (2) are aggregated. Let 𝜎𝜎 be a partition of the column indices {1, … , 𝑃𝑃} into a set of ters𝑆𝑆 𝑘𝑘 , 𝑘𝑘 = 1, … , 𝐾𝐾, such that𝑆𝑆 𝑘𝑘 ∩ 𝑆𝑆 𝑒𝑒 = ∅, ⋃ 𝑆𝑆 𝑘𝑘 = {1, … , 𝑃𝑃} 𝐾𝐾 𝑘𝑘=1

. For the 𝑘𝑘 𝑡𝑡ℎ cluster, let its size as|𝑆𝑆 𝑘𝑘 | = 𝑃𝑃 𝑘𝑘 , so that ∑ 𝑃𝑃 𝑘𝑘 = 𝑃𝑃 𝑘𝑘 . The matrix 𝐴𝐴 𝑘𝑘 is defined to be the submatrix of 𝐴𝐴 consisting of those columns whose indices are in 𝑆𝑆 𝑘𝑘 . Also, the subvectors 𝒄𝒄 𝑘𝑘 and 𝒙𝒙 𝑘𝑘 are defined in the same way. Consider a nonnegative 𝑃𝑃 𝑘𝑘 -vector 𝒈𝒈 𝑘𝑘 = (𝐴𝐴 𝑠𝑠 𝑘𝑘 ), satisfying the following normalizing condition

𝒈𝒈 ∈ 𝐺𝐺 = {𝒈𝒈| � 𝐴𝐴 𝑠𝑠 𝑘𝑘 = 1, 𝐴𝐴 𝑠𝑠 𝑘𝑘 ≥ 0, 𝑘𝑘 = 1, … ,
defines the (column or variable) aggregated problem corresponding to [START_REF] Gilmore | A linear programming approach to the cutting-stock problem[END_REF]. In (3), 𝑿𝑿 is a 𝐾𝐾-vector of aggregated variables. For a given original problem, then, an aggregated problem is determined by the pair (𝜎𝜎, 𝒈𝒈). It is assumed that (𝜎𝜎, 𝒈𝒈) has been chosen so that (3) is feasible. The optimal primal and dual solutions of the aggregated problem are denoted as (𝑿𝑿 * , 𝜋𝜋 * ). Interested readers could consult [START_REF] Rogers | Aggregation and Disaggregation Techniques and Methodology in Optimization[END_REF] for a general framework for aggregation and disaggregation technique and a survey on previous works.

The New Approach

Solving LP (1) takes a long time even using state of the art LP solver CPLEX as depicted in Fig. 2. For a particular instance of LP (1) used to sketch the chart, finding an optimal solution of the LP takes maximum time of 1:08:34, minimum amount of 0:21:34, and in average it takes 36 minutes per year for a complete run for a period of 38 years. Considering that the model should be solved for several consecutive years, at least for 38 years in this project, obtaining the optimal land allocation policy for a long planning horizon becomes quite time consuming. In addition, all the economic data, demand and supply figures along with energy price paths are acquired based on only one scenario out of thousands of scenarios. As a consequence, considering more than one scenario simultaneously for the sake of a stochastic model will result in an even larger problem and will require longer computational time. Therefore, for more comprehensive investigations on optimal land allocations an efficient solution methodology is necessary to solve the large LP problem in a reasonably shorter amount of CPU time.

Fig. 2. CPU time for solving the LP for each year from 2013 to 2050

There are some interesting properties associated with agricultural land allocation that are utilized in our aggregation process. For example, it is quite plausible that nearby land parcels would be allocated to similar activities. The reason for this assumption is that geographically close areas have similar economical and soil properties, and the most profitable land allocation should be almost the same considering accessibility, water resources and costs. Therefore, instead of considering land parcels (cells) of size one square kilometer, we aggregate them into larger chunks of land which are called clusters. Each cluster will act as a decision making unit in the optimisation model, and its economic data are obtained by summing over the data related to all active cells in that cluster.

We start by a lower resolution, say 𝑅𝑅 0 , and divide Australia into a number of large parcels of land. Consider the set of clusters 𝐶𝐶 = {𝐶𝐶 𝑘𝑘 } , where 𝐶𝐶 𝑘𝑘 𝑖𝑖 ∩ 𝐶𝐶 𝑘𝑘 𝑗𝑗 = ∅ for 𝑘𝑘 𝑖𝑖 ≠ 𝑘𝑘 𝑠𝑠 , and ⋃ 𝐶𝐶 𝑘𝑘 = 𝑅𝑅 𝑘𝑘 . Each 𝐶𝐶 𝑘𝑘 is a nonempty subset of 𝑅𝑅. For each cluster, all the economic data including profits and costs of all nine activities and all 24 food products, in addition to other related quantities and capacities are calculated based on related data of active cells in that cluster. Note that in the aggregated LP model the number of land-use constraints is less than in the original problem. The reason is that in the original model there is one land-use constraint for each cell. However, in the aggregated model we are dealing with clusters of cells, and there is one land-use constraint for each cluster. The number of other constraints remains the same as the original problem. An aggregated LP is constructed as explained in section 4 with far few-er variables and constraints, and the optimal primal and dual solutions are obtained by solving the problem using CPLEX.

The plan is to improve the aggregated models consecutively to become as close as possible to the original LP with highest resolution but without solving a large scale problem. Towards this aim, one possible way is to subdivide large clusters and construct smaller ones, and introduce new clusters into the model as new columns. We incorporated reduced costs to separate cells with positive amount inside a cluster from the rest of cells in that cluster. Reduced costs are available at the end of optimisation process. However, these reduced costs correspond to the clusters, and we need to extend this concept to cells inside each cluster.

If, for example, 𝑟𝑟 𝑘𝑘,𝑎𝑎𝑎𝑎 is the reduced cost of the variable that corresponds to agricultural activity in the cluster k, we need to distribute it between the agricultural activity of all cells inside the cluster k so that the total sum is equal to 𝑟𝑟 𝑘𝑘,𝑎𝑎𝑎𝑎 . In other words,

� 𝑟𝑟 𝑘𝑘 𝑐𝑐 ,𝑎𝑎𝑎𝑎 𝑐𝑐∈𝐶𝐶 𝑘𝑘 = 𝑟𝑟 𝑘𝑘,𝑎𝑎𝑎𝑎 .
In this manner, all reduced costs corresponding to all variables related to all activities for each cell in cluster k are calculated. However, before the calculation of reduced costs, we introduce a new concept of semi-reduced costs which is part of the final reduced costs and will be formally defined in the next paragraph. For each cell all semi-reduced costs are calculated for all activities, and the maximum of them is recorded. This quantity represents, to some extent, the most profitable activity for a particular cell. This quantity is also used to distribute the reduced cost of the cluster between all cells proportional to the size of its most profitable activity. The reduced costs are calculated based on semi-reduced costs, and by means of the sign of final reduced costs we can partition each cluster into two sub-clusters containing cells with positive and non-positive reduced costs. The economical quantities for newly constructed clusters are calculated using the information of cells belonging to it, and related information for shrunk clusters are modified considering cells remaining in those clusters. The changes in the old cluster transferred into the LP, and a new column is introduced for each new cluster.

Mathematically, in iteration t, after solving 𝐿𝐿𝑃𝑃 𝑡𝑡 , we obtain optimal primal solution 𝑥𝑥 𝑡𝑡 * and optimal dual solution 𝜋𝜋 𝑡𝑡 * . At this stage, it is easy to calculate the reduced cost of each activity j of each cluster k using 𝑟𝑟 𝑘𝑘 𝑗𝑗 = 𝑃𝑃 𝑘𝑘 𝑗𝑗 -𝜋𝜋 𝑡𝑡 * 𝒑𝒑 𝑘𝑘 𝑗𝑗 . To calculate reduced costs of each cell, perhaps the easiest way is to distribute 𝑟𝑟 𝑘𝑘 𝑗𝑗 uniformly between all cells included in cluster k. However, this approach treats all the cells with equity, and there is not a distinction based on profitability. Another, smarter approach is to distribute 𝑟𝑟 𝑘𝑘 𝑗𝑗 in a way that supports cells with higher profit potentials. In order to calculate reduced cost of each cell in each cluster, first we calculate semi-reduced costs, 𝑟𝑟 𝑐𝑐 𝚥𝚥 �, for the problem in this paper. The semi-reduced cost is the reduced cost without the portion related to the optimal dual value of the corresponding land-use constraint.

In other words, for cell c in cluster k, 𝑟𝑟 𝑐𝑐 𝑗𝑗 = 𝑟𝑟 𝑐𝑐 𝚥𝚥 � -𝛼𝛼 𝑐𝑐 𝜋𝜋 𝑘𝑘 𝑙𝑙 . For a cell c in cluster k, if max 𝑠𝑠 {𝑟𝑟 𝑐𝑐 𝚥𝚥 �} is positive for activity j, the quantity 𝛼𝛼 𝑐𝑐 is defined as

𝛼𝛼 𝑐𝑐 = max 𝑠𝑠 {0, max 𝑠𝑠 {𝑟𝑟 𝑐𝑐 𝚥𝚥 �} ∑ max 𝑠𝑠 {𝑟𝑟 𝑐𝑐 𝚥𝚥 �} 𝑐𝑐 }
After calculation of reduced costs of each cell in each cluster, the cluster is divided into two sub-clusters. One sub-cluster contains all the cells with positive reduced costs, and the other one contains all the cells with non-positive reduced costs. The economic data of each cluster is updated and new clusters are added into the previous aggregated model as new columns. The model is solved again, and the process of dividing clusters is continued until further bifurcation is impossible.

Algorithm 1 Aggregation and Column Generation 1: Make clusters with a low resolution 𝑅𝑅 0 , and calculate economic attributes of each cluster. 2: Build the aggregated LP, and get optimal primal and dual solutions. 3: Calculate reduced costs of each activity for each cell in each cluster, and subdivide each cluster into two clusters based on the sign of reduced costs of containing cells. 4: If no new clusters created (no cells with negative reduced cost) then the solution is optimal. Stop. 5: Calculate economic attributes of new clusters, and update corresponding data in old clusters. 6: Improve the aggregated LP and include new columns for each recently created cluster, and get the optimal primal and dual solutions. 7: If a stopping criterion is not met, go to Step 3.

Computational Results

In this section we present some of the numerical results we obtained through working on application of our algorithm to solve the large scale linear programming problem originated from land use allocation. As seen in Fig. 2 solving a LP problem for each year takes a long amount of time averaging 36 minutes. For a fare comparison, all the numerical experiments are carried out on the same computer. It is worthwhile to mention that the optimisation technique discussed in this paper and all linear programming problems were coded in Python 2.7 language and executed on a computer running a dual core 64-bit Intel(R) Xeon(R) processor at 2.79 GHz with 64 GB RAM.

The linear programming problems were solved using IBM CPLEX 12.5 solver. Table 2 summarizes the output of solving the large scale optimisation problem starting with different resolutions. For our algorithm the stopping criterion is met if we are unable to subdivide the existing clusters further into new, smaller clusters. Resolution 1 is the finest resolution and creates the largest problem. Also, the resolution level of 7569 divides Australia into 417 clusters at the beginning of the computation, and creates the smallest possible LP problem in our experiments.

We recorded starting number of clusters, 𝑃𝑃 𝑠𝑠𝑡𝑡𝑎𝑎𝑟𝑟𝑡𝑡 , final number of clusters, 𝑃𝑃 𝑒𝑒𝑒𝑒𝑒𝑒 , the number of iterations to obtain the final optimal value, and the percentage difference between our optimal value with the optimal value of the original problem for each resolution. The original problem in its finest resolution of 1 has 812,383 clusters, due to considering each cell as a cluster. The algorithm could find the optimal solution in one iteration, as it is impossible to break down clusters. We also record the amount of time taken for CPLEX to find the optimal solution for each resolution in minutes. This column provides us a good measure to compare time performance of our approach. We did not consider the time necessary for the extraction and calculation of parameters in this column, and in the execution time provided in section 5. The data in the table demonstrates some interesting facts about our approach and its capability to deal with large scale linear programming problems. Instead of solving the largest possible problem, the algorithm tries to solve an aggregated problem which is a considerably smaller problem. After each iteration, it tries to subdivide clusters based on their reduced costs. As an example, starting with resolution 7,569 which divides Australia into 417 clusters at the beginning, only in seven iterations and with addition of 67 clusters, the algorithm could reach to the optimal objective function with only 0.49 percent deviation in a small portion of a minute. Comparing this outcome with solving the original problem which needs 53 minutes shows how fast the new algorithm could find the optimal solution. Furthermore, Fig. 3 represents the convergence rate to the optimal solution considering each resolution factor. Combining the information presented in Table 2 and Fig 3 it is apparent that there is a smarter approach of solving the large scale LP in a noticeably shorter amount of time by choosing a lower resolution at the beginning and gradually refining the resolution of clusters in the sequel iterations. With this approach, it is meaningful to build a stochastic optimisation model of the land allocation planning problem, and obtain sensible outcomes in reasonable amount of time and effort.

Conclusion

In this paper we presented a new efficient algorithm to solve large scale linear programming problem originated from optimal land allocation planning. The algorithm combines techniques of column generation, aggregation and disaggregation particularly suitable for this problem. We approximated the large problem with a small LP problem by considering huge parcels of land as a cluster, and then tried to improve the approximation iteratively by subdividing clusters when it is possible. The use of our custom aggregation and disaggregation method allows us to easily solve land use models at a much higher resolution than would otherwise be possible, thus improving the fidelity of the models. In addition, the huge reduction in computational time lays the ground for more complex models. For example an extension to stochastic optimisation incorporating multiple scenarios in a single model is being considered.

Fig. 3 .

 3 Fig. 3. Convergence results for each resolution It is interesting to notice that choosing a lower resolution of 1,600 causes to start with a worse bound than 7,569. However, as shown in the Fig 3 it gets better quick, meaning that less clusters do not necessarily mean worse bounds, and also the smaller clusters are not necessarily all subsets of the larger clusters.Combining the information presented in Table2and Fig 3 it is apparent that there is a smarter approach of solving the large scale LP in a noticeably shorter amount of time by choosing a lower resolution at the beginning and gradually refining the resolution of clusters in the sequel iterations. With this approach, it is meaningful to build a stochastic optimisation model of the land allocation planning problem, and obtain sensible outcomes in reasonable amount of time and effort.

Table 2 .

 2 Results on running the algorithm of different resolutions

	resolution iterations	𝑃𝑃 𝑠𝑠𝑡𝑡𝑎𝑎𝑟𝑟𝑡𝑡	𝑃𝑃 𝑒𝑒𝑒𝑒𝑒𝑒	Difference	t (min)
					(%)
	1	1	812383	812383	0.0	53.598
	9	4	138329	138368	0.16	25.930
	81	5	20402	20482	0.27	0.694
	784	5	2784	2899	0.27	0.020
	1600	7	1608	1742	0.26	0.016
	7569	7	417	484	0.49	0.009


Ralf Denzer
email: ralf.denzer@enviromatics.org

Tobias Kalmes

Udo Gauer

A Software Package for Automated Partitioning of Catchments

Keywords: parallel simulation, hydrological network, graph clustering 1. Precondition

This paper reports about a software package which has been developed to automatically partition hydrological networks (catchments) into clusters of similar size. Such clustering is useful for parallel simulation of catchments on distributed computing systems and is typically based on heuristic graph algorithms.

There have been a few approaches to automatically partition catchments, but literature research indicates that there seems to be no systematic investigation of the usefulness of different graph algorithms for catchment partitioning over a reasonable number of real world data sets. Our study aims at making a step in this direction.

The paper describes the software package, which has been implemented in Java, its pluggable architecture, and initial experiments using the European catchment dataset ECRINS. The paper presents work in progress.

Introduction

Some computational problems of river catchment simulations require large amounts of computing time. It is therefore just natural to aim at parallelizing such computations. At the core of any parallelization, the first question is whether an algorithm can be parallelized and how. As catchment simulations are dynamic flow problems, both spatial and temporal subdivision is an option. A look at the literature suggests that while there have been several proposed solution, a systematic study is lacking. The aim of the work presented in this paper is a long term one: to start with a systematic approach coping with the parallelization of the problem (which means the partitioning of the catchments here) and to end in computational infrastructures for easy deployment of parallelized models.

In order to approach the catchment partitioning problem systematically, we chose to apply a divide and conquer strategy, starting with the most simple problem setup. The experiments reported in this paper are based on the following preconditions:

It is clear that in reality other cases may occur: the graph may contain cycles and the algorithms may be different for different parts of the catchment, but many applications in the literature satisfy the above requirements anyway and binary graphs are common representations for catchments.

For the purpose of this study we have also completely ignored the question when the parallelization of a problem will start to pay off, depending of the compute time per time step per node compared to the amount and latency of data transfer between computing nodes. We were only interested in the first step, the automatic partitioning.

A binary graph is a special form of a directed acyclic graph (DAG). Each graph node has a maximum of two predecessor upstream nodes, and a maximum of one downstream node. The problem of our study is defined as follows:

For a binary graph G, find a clustering

C = { S i , i = 1,n } (1)
into n clusters, where S i are binary graphs and

∪ S i, i=1,n = G (2)
,which means that the clusters taken together represent the graph G. Let |G| be the graph norm (number of nodes) in a binary graph. An optimal clustering C opt is one which satisfies the following equation:

|S i | = |S j |, for all i,j (3)
Equation ( 3) means that for the purpose of parallel computing, it is desired that all clusters have equal size. This reflects precondition 2.

Space and Time Parallelization

A dynamic flow problem can be parallelized in space and time, as long as the flow is represented by a DAG in which downstream nodes only receive input from upstream nodes. There are two principle possibilities to use a distributed computing infrastructure for parallelization: a) a large catchment is split spatially into several catchments, which are computed in parallel, and b) for different time steps or time periods there is a different processor per node. Theoretically one could assign |G| * T (T being the number of time steps) processors to the computing problem, which would however result in unrealistically large hardware needs. In reality it makes sense to allocate one computing node for a clustering of nodes, to compute a time period of reasonable size and then to communicate a partial time series between the clusterings through connecting nodes, in order to balance the trade-off between using multiple processors (positive influence) and communication delays and latencies (negative influence).

The long term interest of this study is to get a better understanding of the different trade-offs, based on a large number of real world datasets, and not based on very few (or only one) dataset or very small (or even toy) datasets. For this purpose it is necessary 1. to spatially cluster many different datasets, including large datasets with many nodes, in order to gain experience how good the spatial partitioning is in realityfor this step a concrete model is not needed, 2. to actually run many model runs based on these clusterings in a distributed environment and measure their performance (which will be quite an effort and will require a lot of automation), or 3. to try to build a theoretical model of the distributed computing system based on execution per time step, average latencies and communication overhead, or 4. to experiment with a dummy model on a real distributed computing system To this end we have started with the first step, and have implemented a software suite for spatial partitioning (or clustering). Along with the core software package, three published algorithms and one new algorithm have been implemented which we have used for the first batch of partitioning experiments.

It is our intention to carry on with the first step with as many datasets as we can get hold of and to implement more algorithms in the future. We also intend to conduct experiments according to step 4. We are also investigating how to move the clustering into a cloud infrastructure in order to provide a dynamic and scalable clustering service.

It should be noted though that this study is carried out with groups of Masters students only, without any external funding sources and support.

Related Work

The literature base on partitioning of catchments is not large but very diverse, approaching the problem from different angles. A 2001 article by M. Grübsch and O. David [START_REF] Grübsch | How to Divide a Catchment to Conquer Its Parallel Processing, an Efficient Algorithm for the Partitioning of Water Catchments[END_REF] gives a good introduction into the problem, its computational complexity and graph properties which may be taken into account when developing partitioning algorithms. They also propose a heuristic algorithm as a solution. In [START_REF] Wang | A common parallel computing framework for modeling hydrological processes of river basins[END_REF], a generalised computational architecture is proposed which allocates computing nodes over a master-slave pattern using load balancing, which is a standard approach in distributed computing. In [START_REF] Yalew | Distributed computation of large scale SWAT models on the Grid[END_REF], a SWAT model is parallelized but it remains unclear how the subcatchments have been partitioned. In [START_REF] Denzer | Parallel simulation of environmental phenomena[END_REF] the authors state that a more generalized approach based on well understood software patterns would help practitioners develop parallel simulation solutions. One commonality of published papers is that whatever is proposed, the examples included in published studies have been small sets of case studies. Our study intends to start a systematic investigation over large sets of cases, starting with the partitioning problem.

The Software Package

The software implementation uses a pluggable architecture. At the core is a class called ClusteringSuite which dynamically manages a clustering run (see fig. 1). It loads clustering algorithms and test criteria (algorithms computing some quality measure of the generated clusters).

A new clustering algorithm or a new test criterion respectively only has to implement a simple interface and will be dynamically loaded based on entries in a property file denoting the algorithms and test criteria to be used. The suite dynamically loads algorithms via a method called addAlgorithm() add test criteria via a method called addTestCriterion(). The package contains two different executable programs, one called Clusterer and one called ExperimentController. While the clusterer will cluster one set of samples given by the user (respectively passed to the program), the experiment controller will generate samples which cluster the same catchment into many different sizes. This is for automation of the clustering experiments.

The Experiment

The Catchment Dataset

For the first study we have used the ECRINS 1 dataset, which can be downloaded from the web site of the European Environment Agency 2 . The dataset contains amongst others a European-wide hydrological network with approximately 1.2 million river segments. The download section of the EEA site contains several datasets. In order to access the hydrological network data, the "EcrRiv" dataset is needed. The dataset is well documented in [START_REF] Eea | EEA Catchments and Rivers Network System[END_REF]. For the experiments we have chosen a variation of catchment size from 1500 nodes to nearly 1460000 nodes (see table 1).

Clustering Algorithms Used

After literature research, 3 algorithms were chosen for the implementation of the first experiments. The choice of algorithms was not easy, as many proposed graph clustering algorithms cope with more general graphs, and the documentation is not always clear and detailed enough to use it as a blueprint for implementation. As a start, we decided to implement three algorithms: K-Means [START_REF] Kanungo | An efficient k-means clustering algorithm: Analysis and implementation[END_REF], RNSC [START_REF] King | Graph Clustering with Restricted Neighbourhood Search[END_REF] and Spectral Clustering [START_REF] Ng | On spectral clustering: analysis and an algorithm[END_REF]. During the course of the investigation, an idea lead to the implementation of a new algorithm called Neighbourhood Clustering. This algorithm was also used in the tests and may be published at a later stage after more experimentation.

Aims

If you want to distribute a catchment simulation over a distributed homogeneous computing infrastructure, it is desirable that all clusters have the same size. Then different computing nodes would not have to wait for each other when they are passing messages between clusters between time steps (or time periods) simulated. In reality this is never achievable with cluster size greater than 1 (one computing node per graph node), because the binary graphs representing real world catchments are highly unbalanced. Real results have clusterings in which clusters are not of equal size. It was our goal to achieve clusterings where the maximum cluster size and the minimum cluster size do differ only by a factor of 2 to 2.5. The rationale behind this aim is that due to the asynchronous nature of operating systems, network communication, latencies and so forth the overall runtime behaviour of a distributed simulation is not predictable anyway, and this aim seemed reasonable, achievable and still practical for distributed simulation. In order to make the results more visible we define the maximum spread as that factor (maximum spread := maxsize / minsize), and the spread curve (figure 2) as a curve which shows the spread for all generated clusters, where the clusters are ordered by size from left to right.

Experiment Setup

In our initial experiments we ran all catchments against several algorithms where it was practical. Spectral clustering, for instance, was not used in all cases as it has considerable runtime requirements. Also clustering the largest catchments is a problem of run-time om standard machines. For each catchment in the experiment, the desired cluster count was chosen, starting with approximately 1/100 of the catchment size and iterating to approximately 1/10 of the catchment size with different step sizes, depending on the size of the catchment (for instance for river Neckar with 1578 nodes, an iteration was performed from desired cluster size 15 to 150, for each integer in between. These initial experiments were aimed at getting a first idea what might be achievable with the algorithms chosen.

Discussion

Our initial findings are that not all results are in the order of magnitude which we had hoped. Some clusterings have a very large variation in number of nodes. The best spreads are in the order of magnitude of 2.5, but particularly the very large catchments seem to be difficult to cluster evenly. Some algorithms produce stray nodes, clusters of very small size -in some cases many of them. It is not clear yet whether those can be easily re-connected to other clusters (or to each other) in post-processing without unbalancing the quality of the solution again.

We have also not quality-checked the implemented algorithms over a large number of experiments, and a real problem is that there are no datasets in the literature to which we could compare our results to.

Also the choice of starting points influences the clustering, and some algorithms do have parameters which influence their behaviour. We have not yet done any variation over the automatically chosen starting points and over these parameters.

The computation time for the clusterings is another problem. It is just not practical to run many clustering runs over reasonable size clusters on standard machines, as some algorithms (particularly Spectral Clustering) just need too long to compute one clustering. Therefore particularly the large catchment have not been run over all variations of algorithms and cluster sizes which we had initially planned.

Therefore, in another project course in 2014/2015, we plan to extend our solution to a cloud infrastructure, in which we aim at running many more variations of experiments in order to find better solutions.

Fig. 1 .

 1 Fig. 1. Parts of the clustering suite interfaceThe current implementation provides the following test criteria as built-in elements:• AverageSize gives the average size of the clusters • AboveAverageSize counts how many clusters are above average • BelowAverageSize counts how many clusters are below average • MedianSize gives the median size of the clusters • AverageWeight gives the average weight of the clusters • MedianWeight gives the median weight of the clusters • DeltaMedianAverageSize gives the delta of average size and median size • Correctness determines whether the algorithm is correct at all

Fig. 2 .

 2 Fig. 2. Real vs. ideal spread curve (Neckar, 1500+ nodes, NH algorithm, 15 clusters)

Table 1 .

 1 The river id is the column river_id in table c_tr of ECRINS and the end node is the draining node (column nodid in table c_node of ECRINS). This data was run against several algorithms. Catchments used in initial experiments

	River	River id	End node	Nodes
	Volga	Z_C0000603	Y000601274	146585
	Danube	Z_C0000897	Y000828008	114318
	Kama	Z_C0000605	Y000312116	64024
	Dnieper	Z_C0000089	Y000728340	42962
	Don	Z_C0000631	Y000640649	32398
	Rhine	Z_C0000823	Y000617080	25270
	Po	Z_C0000025	Y000962481	23074
	Ural	Z_C0000588	Y000537580	21391
	Ebro	Z_C0000053	Y001234399	15648
	Douro	Z_C0000049	Y001129129	12459
	Drave	Z_C0000933	Y000912526	11966
	Loire	Z_C0000067	Y000774881	9815
	Adige	Z_C0000088	Y000952540	7608
	Isere	Z_C0000750	Y000956258	6954
	Adda	Z_C0000020	Y000955093	3385
	Moselle	Z_C0000899	Y000660011	3054
	Neckar	Z_C0000896	Y000693135	1578
	1 http://www.eea.europa.eu/data-and-maps/data/european-catchments-and-rivers-network
	2 http://www.eea.europa.eu/

Acknowledgements

This experiment was carried out as part of a Master level project course in 2013/2014, with a group of 6 students, of which 2 agreed to co-author this paper. The resulting software was cleaned and refactored by the main author after the course had produced the first results. It is our intention to continue these investigations and we welcome collaboration with external partners.


Rodolfo García-Flores

Magnus Moglia
email: magnus.moglia@csiro.au

David Marlow
email: david.marlow@wiseranalysis.com

Rodolfo Garcia

A Framework for Optimal Assessment of Planning Investments in Urban Water Systems

Keywords: urban water planning, infrastructure planning, decision making for water networks

The functionality expected by governments and citizens from urban water management systems (UWMS) has evolved in time from delivering basic services to enabling complex issues such as healthy ecosystems, environmental sustainability, and economic growth. Alongside these changing expectations, the pressure on policy makers to fulfill disparate performance metrics with diminishing resources is increasing. In this paper, we introduce the Water Assets and Infrastructure Network Decision Support (WAND) framework for modelling and planning urban water systems with the aim of assisting medium to long-term investment decisions. The framework combines economic performance measures with liveabililty and sustainability, and uses optimisation as a crucial component. To demonstrate its feasibility, we present a prototype environmental software system built upon the principles of legibility, adherence to engineering conventions, and extensive unit testing. We use the prototype system to analyse a model of a hypothetical urban water system formed by two coupled water networks, one for freshwater and one for stormwater collection, and which handles six different commodities. Our results suggest to planners the optimal combination of planning investments while considering capacities, service levels and network operating conditions.

Introduction

The Water Services Association of Australia (WSAA, 2013) estimate that the urban water sector delivers services to over 20 million Australians in more than 9 million connected properties. It manages over AUD $120 billion in assets and between 2006-07 and 2011-12 capital expenditure was estimated at over $33 billion. Operational expenditure in 2011-12 was estimated at $7.2 billion. The urban water sector is a significant part of the Australian economy. More than that, urban water infrastructure represents a significant capital stock and future generations will inherit the cumu-lative impact of our on-going investment decisions [START_REF] Nolan | Funding Water and Wastewater Projects: Growing Communities Seek Strategic Alternatives[END_REF][START_REF] Marlow | A pathway to a more sustainable water sector: sustainability-based asset management[END_REF]. The core urban water services involve the management of three water streams: water supply, sewage (and used water), and stormwater. Traditionally, each of these streams has been managed separately, but pressures on these systems have, however, grown over time and there is an increasing clash between the demand for and limits of water resources. This clash results in ecological, economic and cultural 'strains' (Vlachos and Braga, 2001), which have led various authors to suggest that the current model of water service provision is not sustainable, especially when considered in light of future uncertainties such as climate change and population growth [START_REF] Ashley | Sustainable decision making for the UK water industry[END_REF][START_REF] Milly | Climate change: Stationarity is dead: Whither water management?[END_REF]; [START_REF] Wong | Transitioning to Water Sensitive Cities: Ensuring Resilience through a new Hydro-Social Contract[END_REF].

From a practical perspective, innovations are being introduced into systems dominated by traditional centralised solutions to produce hybrid systems that combine decentralized green infrastructure systems, such as storm water harvesting, and traditional capital-intensive infrastructure. Balancing the costs, benefits and risks of both approaches requires enhanced capacity to predict the future and consider multiple tradeoffs. This to a large extent becomes an optimisation problem, which involves effective management of infrastructure assets. Further discussion of such management strategies is provided by [START_REF] Maheepala | Integrated Urban Water Management Planning Manual[END_REF] or [START_REF] Marlow | Towards sustainable urban water management: A critical reassessment[END_REF].

Regarding optimal design of water networks, Daganzo (2012) considered public infrastructure systems where the policy-making body chooses prices and system design by maximising society's welfare. This combined problem was decomposed into sub-problems, and optimal operation was induced by structuring payments adequately. [START_REF] Bieupoude | Optimization of drinking water distribution networks: Computer-based methods and constructal design[END_REF] used a geometric optimisation method known as "constructal" design to optimize water distribution network architectures subject to operational water quality constraints. [START_REF] Marinoni | Directing urban development to the right places: Assessing the impact of urban development in an estuarine environment[END_REF] developed a method based on mathematical programming to determine the impact of development on a waterway that receives treated water based on water quality criteria, and applied it to Greater Sydney Metropolitan Area.

There exist environmental software systems that inform investment in UWMS, including a small number of commercial products which to inform optimal investments in urban water systems. However, none of these consider all the metrics simultaneously. Among the most widely used are PIONEER1 [START_REF] Lumbers | Optimal asset management planning: advances in water mains and sewers analysis within a new modelling environment[END_REF], and Wil-CO 2 [START_REF] Burns | Hydrologic shortcomings of conventional urban stormwater management and opportunities for reform[END_REF], which are suites of models built as environments for planning water infrastructure. Both WiLCO and PIONEER have been developed in the UK institutional and regulatory context in mind, though their parent companies claim that they can be extended to the requirements of other sectors. The modelling undertaken within WILCO and PIONEER focuses on water and wastewater assets, which are considered in isolation from each other, and trade-offs are thus analysed within each service area, but not across service areas. Furthermore, neither WiLCO nor PIONEER were designed to represent the urban water system as a whole, so they do not consider wastewater, drinking water and stormwater in an integrated manner.

In this paper, we introduce the Water Assets and Infrastructure Network Decision Support (WAND) framework, which has the aim of assisting medium-to long-term investment decisions, and relies on optimisation as a crucial component. The proofof-concept solves a mathematical model of a hypothetical UWMS. Our results suggest to planners the optimal combination of investments, while considering capacities, service levels and operating conditions for sets of equipment and facilities.

2

The WAND Framework

The WAND framework (Figure 1) is based on the fact that baseline performance as reflected in the system of representative metrics (related to liveability, sustainability and productivity) will vary over time due to natural changes, including the impact of asset deterioration and other natural factors such as of climate change. Performance will also be influenced by anthropogenic factors such as changes in customer requirements/expectations, regulatory requirements, population growth or macroeconomic development. The impact of these changes on system performance can only be countered by making investments in capital and non-capital interventions; informing these is the purpose of the core engine of WAND, which relies on optimisation to reduce the gap between actual and target performance. To feed the optimiser, information is required on the policies and strategies to implement and the effectiveness of the associated intervention options. The optimal investment decisions, selected from a feasible option space, alter the system as part of a feedback loop that operates at strategic, tactic and operational time scales. This helps form a self-similar model, in the sense that it is as representative to consider specific service areas as it is to consider the complete urban water system. In fact, in theory the framework could be expanded to include all urban systems at the three key decision scales and provide systemic insights into the impact of investments across a whole range of service provision areas.

The WAND prototype software system introduced in this paper implements the tasks of the core engine, and its purpose is to assess the tradeoffs on livability, sustainability and productivity associated with different portfolios of investment within a single framework. The prototype handles investments by distinguishing between asset operation, replacement, maintenance and construction.

To illustrate the purpose of the WAND environmental software system, we address on the problem of urban water runoffs. As cities grow the landscape changes and there is an increase in the area of land that is covered by impervious surfaces. This means that, when rain falls on the city, much of the water that would otherwise have penetrated into the soil now gathers as urban runoff. This changes the flow patterns in rivers and creeks, with further impacts on ecosystems and erosion of river banks, as well as increasing the risk of urban flooding with consequent societal costs. Furthermore, the runoff in cities goes directly to receiving waters without treatment [START_REF] Burns | Hydrologic shortcomings of conventional urban stormwater management and opportunities for reform[END_REF], becoming a significant source of pollution. WAND addresses runoffs by assisting investment decisions in:

• Capture and use of stormwater; i.e. stormwater or rainwater harvesting [START_REF] Sharma | Environmental and Economic Assessment of Urban Water Services for a Greenfield Development[END_REF]. This includes the treatment of urban runoff as per local regulations. • Water Sensitive Urban Design features [START_REF] Wong | An Overview of Water Sensitive Urban Design Practices in Australia[END_REF], which are helpful in order to achieve a more natural water cycle and provide additional treatment of urban runoff, such as stormwater bioretention systems and constructed wetlands. These are also known as green infrastructure.

Computational System

The optimiser in the WAND framework was implemented in Clojure3 1.5 and CPLEX 4 12.4 with an Excel interface, as the UML 2 component diagram [START_REF] Jacobson | The Unified Software Develop-ment Process[END_REF] in Figure 2 shows. Clojure is a compiled dynamic programming language of the LISP family built for the Java Virtual Machine (JVM) that has features that make it an attractive alternative to code the solution of problems that depend crucially on a correct mathematical representation. First, it has all the flexibility and power of LISP, but without the burden of backward compatibility with historical versions (Halloway, 2009). Second, it connects seamlessly to the mainstream infrastructure provided by Java by embracing the JVM. It does not attempt to build wrappers around Java, but tries to access it idiomatically, enabling the programmer to leverage on his past programming experience and on the huge amount of legacy code available for Java. Additional advantages include good concurrency support, a small core, laziness and conditions for error handling. Our WAND implementation makes extensive use of the symbolic features of LISP, for example for indexing of variables, making the order of the indexes irrelevant. This is extremely useful because the client's expert knowledge is often based in the particular situations of individual water management resources, and using names is easier to understand than referring to elements by looking up tables of indexes.

Fig. 2. UML 2 component diagram of the WAND optimiser

We developed the system with clarity as the core requirement. We aimed at keeping the problem as close to its mathematical representation as possible. To achieve this goal, we incorporated the following development principles:

1. Literate approach. When coded in literate style, a program should explain a problem and its solution to a human reader because its concepts have been introduced in an order that is best for human understanding [START_REF] Knuth | Literate programming[END_REF]. This implies that software writing has more to do with the craft of writing, including its gradual process of improvement and refinement of ideas, than with merely laying down instructions for a machine to perform. The project was designed with namespace WAND::weave to process the documentation by weaving as introduced by [START_REF] Knuth | Literate programming[END_REF]. All documentation in Clojure is stored as metadata, which in our system is retrieved by the functions in this namespace and "woven" (i.e., formatted) into L A T E X source. The documentation in the woven document is divided in sections, each of which corresponds to a namespace in the optimiser. Each section includes a description of the namespace's functions and variables. L A T E X enables the inclusion of formulas or typesetting commands as part of the documentation. 2. Unit test design. Unit testing takes advantage of the structure that Clojure gives to metadata. The :test field of the metadata enables the storage of a test function as part of a variable. By retrieving the names of all functions (i.e., data members that have argument lists) in a namespace, it is possible to generate lists of functions to be tested automatically, greatly facilitating unit testing.

An extensive and automated unit testing namespace was implemented in WAND::test-support. 3. Dimensional analysis. In our experience, the lack of support for dimensional analysis has been the source of multiple problems in the past. We decided to take advantage of existing code from the JScience5 library. Re-using Java code is encouraged in languages developed for the Java Virtual Machine, and is a particularly simple task in Clojure. Because dimensional analysis is one of the most basic engineering tools, the selected framework should be easy to reuse (namespace WAND::dimensions). Namespace WAND::preprocess translates Excel input data into a data structure that can be handled by the optimiser. This is used by WAND::problem to produce the mathematical formulation, using the support functions stored in WAND::create-lp. Namespaces WAND::log and WAND::display keep a log of calculation progress and format the input and output data as KML maps, respectively. For the external components, the Excel interface captures input data regarding existing and proposed infrastructure, demand, intake levels and capacities, and retrieves the results in spreadsheet format. It also has the capacity to display input and solutions as maps in KML format, by invoking Google Earth 6 . CPLEX is a commercial solver for mixed integer linear programming problems.

Although there is great variety of tools available for development of environmental software systems [START_REF] Argent | Land and water resources model integra-tion: software engineering and beyond[END_REF], we tried to keep the number of tools to a minimum to enhance clarity. The structure of the WAND optimiser in Figure 1 is the product of assimilating previous and successful experiences in coding similar optimisation problems (Higgins et al.; 2013, García-Flores et al., 2014) using the above-mentioned coding principles. The main advantages we perceive from the approach to software development adopted for WAND include shorter length of source code, legibility, reliability and facility to produce documentation. Most importantly, because the code is close to the notation of the mathematical formulation, we also expect WAND's optimisation module to be refined further and reused in future projects.

Example -Melbourne Urban Water System

This example illustrates the implementation of the framework to the Melbourne water system (Figure 3), which consists of two coupled water distribution networks.

• The fresh water network, which also incorporates wastewater streams, is shown in solid arrows and transports fresh water (that is, water for indoor and outdoor use) from reservoirs (sources marked with s f in Figure 3) to the clients in Melbourne's suburbs (marked r). After use, the waste water collected is sent to treatment plants, marked u f , and part of it is used again, (this is shown as double arrows to indicate separate pipelines) and part is disposed of in Port Philip Bay (shown as a gray block arrow). The "return" direction of the double arrows represents what is known as "purple pipes", which is a separate network of limited reach. Most areas do not have purple pipes going to them. • The stormwater collection network, shown in dashed arrows, collects rain water in catchments s m in different points in the city. After some treatment in sites marked as u m (wetlands, for example), part of it is used, depending on the setup, as indoor or outdoor water (or indoor only or outdoor only), and part of it ends in Port Philip Bay. One of the objectives of the treatment facilites u f and u m is to reduce the amount of nutrients thrown into the bay. Irrigation areas can be considered sinks, and green infrastructure is modelled as a subset of treatment facilities u m , which comprise wetlands, vegetated swales and biofiltration. The success of these facilities depends on their ability to store, retain or detain spill-overs, which means that they should also be modelled as reservoirs. Green infrastructure has a positive impact on stormwater flows in different areas. In the model, sewer overflows need to be penalised, whereas stormwater overflows need to be managed.

Fig. 3. The Greater Melbourne Water System

The WAND model considers six commodities and their transformations:

• Indoor water is mostly potable water, but some of the indoor water use may be covered by non-potable sources, i.e., not treated to class A standard. • Outdoor water is mostly irrigation water. It travels from treatment facilities u m to clients r, and from treatment plants u f back to clients r, in both cases to be used as outdoor water. In the path from u f to r, outdoor water is really treated waste water. The final destination of outdoor water is the client nodes. • Fresh water is an aggregate of indoor and outdoor water. It travels from sources of fresh water s f to households r. At the households, the water becomes indoor or outdoor, depending on its use. • Stormwater is sent from stormwater catchments s m to treatment facilities u m .

Stormwater can only be used for outdoor purposes. • Waste water is sent from clients r to water treatment plants u f . • Discharge water is the water returned to the environment into the bay r b . The problem is formulated as a multi-period mixed-integer linear program whose aim is finding the optimal topology and capacities of pipelines connecting the nodes of the water network. This is done by minimising the total costs (i.e., the objective function), expressed as the sum of operation costs (including spillage costs and the costs of exceeding the control limits of sites that act as reservoir), replacement, maintenance and construction costs of all infrastructure assets.

Along the time horizon, the solution of the problem provides information about facility and pipeline segment replacement, upgrade and new construction. The constraints of the problem describe the conservation of flow at facilities and the conversion between water commodities. Some of the facilities are modeled as reservoirs, such as the waste water treatment facilities and green infrastructure (e.g., wetlands). These facilities need constraints with an upper capacity limit, which, if violated, (that is, if the facility overflows) are penalized in the objective function. Sources and stormwater catchments are facilities that provide an influx of fresh and stormwater, respectively, into the system, based on experts' forecasts and which also act as reservoirs, but are not penalised for spillage.

Results and Discussion

The planning horizon considered was from 2014 to 2019, in annual periods, and the tested scenarios were:

1. The creation of a new suburb in 2017 in node 009 with the consequent increase in fresh water demand, together with a linear increase in demand of node 011 from 481 ML of indoor use water per year in 2014 to 1000 ML in 2019 in node 011. All pipeline segments start with a diameter of 150 mm, except for the segments connecting the network to node 009, which do not exist at the beginning of the time horizon. 2. A "business as usual" scenario where no new suburb is projected and no increase in demand for node 011 is foreseen.

The demands of the clients to test in the scenario are shown in Table 1; the assumed available technologies and their capacities, which represent the investment options, are listed in Table 2. The calculation parameters are listed in Table 3. We obtained all the following results in a 64-bit Intel Xeon CPU with 2 processors of eight cores (2.27 GHz) each and 48 GB of RAM. The linear program of the preprocessing stage problem has 22,542 real variables, of which 14,880 are binary and 7,762 are real, and 50,421 constraints. A typical run takes around two hours and forty minutes to complete. Tables 4 and5 show the suggested investments under scenario 1. These are, first, the construction of a fresh water pipeline and on-site infrastructure for node 009 to fulfill the new demand, as expected. Second, the increase in capacity of the pipelines connected to the waste water treatment plant in node 012. These also occur under scenario 2, as the initial pipe diameter assumed for the flows to node 012 is insufficient. Third, we can see the increase in capacity of the pipelines that connect node 011 to node 013, due to the increase in demand in node 011. Finally, we note in Table 5 that the optimiser suggests a "downgrade" of the fresh water source in node 002. This is probably not an action that a policy maker would take, but in this case we obtain this result because the cost of downgrading is not high enough. The decision maker can easily change the value of this parameter. Finally, scenario 2 turns out to be 11.3% more costly than "business as usual". The expectations from urban water management systems have grown to include not only the provision of basic services, but also enabling healthy ecosystems, environmental sustainability, and economic growth. To achieve these goals, we have introduced the Water Assets and Infrastructure Network Decision Support framework. WAND uses optimisation and is built using the principles of legibility, adherence to engineering conventions, and extensive unit testing. To demonstrate its feasibility, we applied a prototype software system to a model of an urban water network that handles six different water commodities. The results suggest to planners the optimal combination of planning investments, while considering capacities, service levels and network operating conditions. The main advantages we perceive from using standard engineering conventions (such as dimensional analysis) and the coding principles of literate programming and unit test design are, first, a shorter length of source code when compared to similar C++ developments; second, it is easier to produce documentation; and third, increased legibility and reliability. The main strength of the WAND optimiser is that it attempts to keep source code as close to the mathematical representation of the problem as possible by taking advantage of LISP's flexibility. This facilitates debugging and enables new members of the development team to "hit the ground running". On the downside, the disadvantages of our approach are the slow execution times when compared to C++ code and the unavoidable learning curve if the analyst is not familiar with LISP syntax. Nevertheless, we expect WAND's optimisation module to continue being refined and reused in future projects, although concepts such as deterioration are yet to be incorporated.

Fig. 1 .

 1 Fig. 1. The Water Assets and Infrastructure Network Decision Support (WAND) framework.

Table 1 .

 1 Demands of clients in the Greater Melbourne case study in ML/year

	Node (clients) Indoor water demand	Outdoor water demand
	Node 007	1326.151	241.092
	Node 008	648.180	90.411
	Node 009	165.758 (after 2017)	23.121 (after 2017)
	Node 010	117.47	16.385
	Node 011	481.356 (increasing to 1000 by 2019)	67.142

Table 2 .

 2 Available technologies to install in nodes and pipeline segments

	Infrastructure type	Available technologies
	Nodes	Site 200 ML, Site 500 ML, Site 1000 ML, Site 1500 ML.
	Pipeline segments	Pipe 150 mm, Pipe 250 mm, Pipe 350 mm

Table 3 .

 3 Penalisation parameters to enforce supply and operational levels in reservoirs

	Penalty for failing	Penalty for fail-	Penalty for violating	Penalty for overflow-
	to deliver indoor	ing to deliver	control limit in storm-	ing stormwater col-
	water	outdoor water	water collection point	lection point
	$10.0M /ML	$10.0M /ML	AUD $1000 /ML	AUD $0.1M /ML

Table 4 .

 4 Results for the upgrade and construction of new pipeline infrastructure

	Origin	Destination	Commodity	Year Old infrastruc-	New infra-
	node	node		ture (mm)	structure (mm)
	012	007	Fresh water	2015 150	350
	012	008	Fresh water	2015 150	350
	013	011	Fresh water	2015 150	350
	007	012	Waste water 2015 150	350
	008	012	Waste water 2015 150	350
	011	013	Waste water 2015 150	350
	004	009	Fresh water	2017 None	350

Table 5 .

 5 Results for the upgrade and construction of new pipeline infrastructure

	Node	Node type	Commodity	Year	Old infrastruc-	New infrastruc-
					ture (ML)	ture (ML)
	002	Source	Fresh water	2016	1500	500
	002	Source	Fresh water	2017	500	200
	004	Junction	Fresh water	2017	None	200
	009	Client	Fresh water	2016	None	200
	6	Summary, Conclusions and Future Applications

http://www.tynemarch.co.uk/products/pioneer/index.shtml, accessed on February

28, 2014. 2 http://www.seamsltd.com, accessed on February 28, 2014.

http://clojure.org/, accessed on the 8 of April of 2013.

http://www-01.ibm.com/software/integration/optimization/cplex-optimizer/, accessed on the 8 of April 2013.

http://jscience.org/, accessed on the 30 of October 2014.

https://www.google.com/earth/, accessed on the 30 of October 2014.


Marshall Ward
email: marshall.ward@anu.edu.au

Yuanyuan Zhang
email: yuanyuan.zhang@au.fujitsu.com

Scalability of Global 0.25°Ocean Simulations Using MOM

Keywords: ocean modelling, performance profiling, high performance computing, parallel computing

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Current global climate simulations typically rely on coarse ocean models of approximately 1° resolution, but there is growing demand for a greater resolution of the ocean eddy fields, with corresponding resolutions on the order of 0.25° or 0.1° [START_REF] Taylor | An Overview of CMIP5 and the Experiment Design[END_REF]. The very strong stratification of the ocean causes its most turbulent currents to emerge at these smaller resolutions [START_REF] Chelton | Geographical variability of the first baroclinic Rossby radius of deformation[END_REF], which are absent from coarse-resolution climate models. Such turbulent processes are expected to have a major role in the maintenance of the ocean's strongest currents, and in the vertical mixing and stratification of the ocean [START_REF] Farneti | The Role of Mesoscale Eddies in the Rectification of the Southern Ocean Response to Climate Change[END_REF].

Along with the many scientific challenges of high-resolution climate modelling, an additional impediment in the adoption of greater resolutions in the ocean is the significant computational cost. A typical climate simulation requires decades, if not centuries, of simulation time to achieve a scientifically valuable result, and can require hundreds of such runs to perfect the tuning of their numerous parameters [START_REF] Taylor | An Overview of CMIP5 and the Experiment Design[END_REF]. In such an environment, one needs the capacity to simulate many years per day. To achieve such results will require proven scalability into hundreds, if not thousands, of CPUs.

In this paper, we assess the ability to run high-resolution simulations on computing platforms of Australia's National Computational Infrastructure (NCI). We focus on scalability up to 1920 CPUs on the Raijin and Fujin computing platforms. We consider a selection of configurations on each platform and assess their efficiencies for use in future scientific research.

Methods

Numerical Model

The experiment used in this study is the global ocean-sea ice model of [START_REF] Spence | Rapid subsurface warming and circulation changes of Antarctic coastal waters by poleward shifting winds[END_REF], which was based on the Geophysical Fluid Dynamics Laboratory (GFDL) CM2.5 model [START_REF] Delworth | Simulated climate and climate change in the GFDL CM2.5 highresolution coupled climate model[END_REF]. The numerical submodels are the Modular Ocean Model (MOM) and the Sea Ice Simulator (SIS), built from the Flexible Modeling System (FMS) framework, which also includes the ocean-ice coupler [START_REF] Griffies | Elements of the Modular Ocean Model (MOM)[END_REF]. The simulations presented here use the MOM 5.1 source code release, which includes the SIS and FMS components.

The numerical grid of the ocean and sea ice models is a curvilinear grid with a nominal resolution of 0.25° and contains 1440 × 1080 horizontal grid points. The ocean and sea ice models use 50 and 6 vertical levels, respectively.

Parallelisation is achieved by decomposing the horizontal grid into tiles of equal size, with additional halo grid points surrounding each tile. For each timestep, the halos of adjacent and diagonally adjacent tiles are updated, requiring 8 messages per field.

Standard experiments are run for 31 days, using 1488 timesteps of 1800 second resolution. Each ocean timestep requires 72 additional sea ice timesteps. A longer simulation time was chosen to represent a typical integration time for scientific analysis, and to reduce any statistical variability of the observed walltimes. Simulations using 960 CPUs were run three times to confirm the reproducibility of the results, and the median walltimes was used to select the run in this study. For all other runs, a single simulation was used for each configuration.

Diagnostic output consists of 4 three-dimensional fields, which are saved to disk after every 5 days. This output rate was chosen to reflect a comparable rate of scientific value while also not overshadowing the model calculations.

Raijin

Raijin is the principal supercomputer of the NCI National Facility. It is a Fujitsu PRIMERGY cluster comprised of 3592 computing nodes, with each containing two 2.6 GHz 8-core Intel Xeon Sandy Bridge (E5-2670) CPUs, with a total core count of 57472. Turbo boost is enabled for these runs, increasing the maximum clock speed of 16 cores to 3.0 GHz. Operational nodes have approximately 32 GiB of memory, with more memory available on selected nodes.

The cluster network uses InfiniBand (IB) FDR-14 interconnects, which provide peak (4x aggregate) transfer speeds up to 56 Gb/s between nodes. There are 18 nodes connected to each top-of-rack (ToR) IB switch. The 72 nodes in each rack use 4 ToR switches. Each ToR switch connects to a top-level director switch, with 3 connections per ToR to each director. Because of the high usage rate of Raijin, it is prohibitive to arrange jobs within the network, and we assume that all experiments in this study depend on the director switches.

The operating system platform is CentOS 6.5, which uses the Linux 2.6.32 kernel. Jobs are managed using the PBSPro scheduler, and files are kept under a Lustre 2.5 filesystem.

We use the Intel Compiler Suite 14.3.172 for model compilation. MPI communication uses Open MPI 1.8.2 and data IO uses NetCDF 4.3.0, which was built using HDF5 1.8.10. The model is built with the -O2 and -xhost optimisation flags.

We consider three experiment configurations on Raijin: a standard configuration using 16 model processes per node (herein "default"), a similar configuration using hyperthreaded cores, and a third configuration running 12 model processes per node (herein "12 PPN"). For the 12 PPN configuration, turbo-boosted clock speeds are increased to 3.1 GHz.

Fujin

Fujin is NCI's Fujitsu FX10 cluster consisting of 96 compute nodes, each containing a 1.848 GHz 16-core SPARC64 IXfx CPU. For our investigation, numerical experiments were limited to a maximum of 84 nodes, or 1344 cores.

The interconnect is a streamlined version of the K-Computer's Tofu interconnect system, and consists of 8 serially connected Tofu units. Each unit contains 12 compute nodes arranged as a 3D torus of shape 2× 3 × 2. Nodes are connected by 10 bidirectional links with 5 Gb/s transfer rates, yielding a peak transfer rate of 100 Gb/s between nodes. Tofu interconnects also include the Tofu barrier hardware module for optimised MPI reductions and broadcasts, but we do not use it in this study.

When compared to Raijin, the lower clock speed of Fujin's CPUs will necessarily result in lower optimal performance for computationally-bound software. However, the fixed Tofu interconnect offers the potential for optimized node layout and efficient scalability for higher CPU counts. The FX10 operating system is Fujitsu XTC OS 1.3.6, with a Linux 2.6.25 kernel. Jobs are managed by the scheduler of the Parallelnavi suite. The filesystem is the Fujitsu Exabyte File System (FEFS) 1.3.1, which is based on the Lustre 1.8.5 filesystem.

The model is compiled using the Fujitsu compiler suite 1.2.1 and includes Fujitsu's MPI library, which is based on Open MPI 1.4.3. Data IO uses NetCDF 4.3.2, built with HDF5 1.8.12 and NetCDF-Fortran 4.2. The model is compiled using the -Kfast, -Kocl, and -O2 optimization flags.

We consider two configurations on Fujin, denoted here as snake and ladder layouts. Node placement for a single Tofu unit for each layout is illustrated in Figure 1. The snake layout guarantees that all communications is exclusively between neigh-bour nodes, but also requires that each node span an entire latitude line. The ladder layout divides each latitude line across two nodes, allowing for smaller tile widths, but it also requires diagonally adjacent tiles to send data through non-neighbour nodes.

Timing

On Raijin, we use the IPM 2.0.2 profiler to construct basic runtime statistics. Model runtime and MPI usage times are based on IPM output logs. On Fujin, runtime and MPI usage times are from the fipp profiler output logs. Submodel runtimes are provided by the internal timers of the MOM ocean model, which are based on the system_clock subroutine of the Intel Fortran compiler and the system clock of the Linux kernel.

Results

Platform Configuration

Figure 2a shows the total runtime for a 480-CPU experiment across different hardware platforms and configurations. The configuration choice on Raijin had a major impact on performance, with observed speedups of 28% and 60% for hyperthreaded and 12 PPN experiments, respectively. Process layout on Fujin had an observable but less notable impact, with the snake layout achieving a 15% speedup in comparison to the ladder layout. Experiment runtimes on Raijin are faster than on Fujin, with the fastest runs on Raijin outperforming the Fujin jobs by more than a factor of two. Much of this difference in performance can be attributed to the clock speeds of the two platforms. However, the clock speed ratio of the Raijin and Fujin CPUs is only 1.62, and cannot alone account for the differences in performance. where f is the CPU clock speed, p is the mean fraction of MPI instructions across all ranks, and τ o is the ocean submodel runtime. We focus on the ocean submodel because it is strongly dominated by floating point arithmetic, and is shown in the next section to demonstrate high scalability up to 1920 CPUs.

N c is an imperfect estimate of cycle count, since total communication time may not accurately represent communication time in the ocean model. There may also be addi-tional communication costs outside of the MPI library. But it can provide an approximation of the number of computational cycles, as well as the relative efficiency of vector arithmetic on each platform.

After correcting for clock speed and focusing on the ocean model, our estimate indicates that the performance of the two platforms is comparable, with SPARC operations requiring only 15% more cycles than the x86 operations for its simulation.

Scaling

The general performance and scalability of each experiment across different CPUs is shown in Figure 3. The general trends shown in the previous section are also observed in Figure 3a, with simulations on Raijin generally outperforming those on Fujin. Hyperthreaded and 12 PPN jobs continue to outperform the default runs on Raijin, and snake layouts outperform ladder layouts on Fujin. The most notable observation is the dramatic loss of performance of the default experiments on Raijin, where scalability begins to diminish at 480 CPUs and walltimes begin to increase after 960 CPUs. Scalability is restored by either introducing hyperthreading or by using fewer cores, indicating the presence of a computational bottleneck related to the number of processes per node.

The relative efficiency of each configuration is shown in Figure 3b, which is computed as ε= rτ 𝑟𝑟 nτ 𝑛𝑛 where n is the number of CPUs, τ n is the runtime using n CPUs, and r is the reference experiment, which is 240 CPUs for all configurations except the ladder Fujin runs, which use 255 CPUs.

Figure 3b shows that, with the exception of the default Raijin setup, all runs demonstrate scalable performance up to 960 CPUs and maintain an efficiency greater than 80%. Experiments on Raijin scale slightly better than on Fujin, although the results are comparable over all platforms and configurations. Beyond 960 CPUs, all experiments begin to drop substantially and no longer provide a reasonable level of performance.

Model performance can be further clarified by investigating the scalability of the respective submodels, as shown in Figure 4. Figure 4a shows that, despite a reduction of general scalability around 960 CPUs, the ocean submodel continues to scale efficiently to 1920 CPUs. However, the sea ice and coupling model subcomponents, shown in Figures 4b and4c, grow progressively worse in performance as the number of CPUs increase, indicating major bottlenecks in the sea ice model and submodel flux exchange components.

There is some evidence of improved efficiency in the ocean model, indicating that the communication bottlenecks in the sea ice and coupler models may also be freeing additional resources for the ocean model. Scaling efficiency on Fujin appears to be better than Raijin for the ocean model, although it is notably lower in the other model components.

Discussion

Configuration

Scalability on Raijin beyond 480 CPUs was only possible after either reducing the number of processes per node, or enabling hyperthreading on the core. In both situations, we reduced the number of processes per computational core. While the underlying cause is not known, one point of consideration is use of subthreads within Open MPI. Since each model instance create four MPI subthreads, there is a competition for resources between these processes on each core. By allowing the kernel to shift some of these jobs to other computational cores, we may be relieving any bottlenecks related to shared memory management or context switching. But further investigations are required to confirm this explanation, and to determine if it is a consequence of x86 architecture, kernel scheduling, or Open MPI implementation.

Scalability on Fujin was not affected by processor layout, although the performance of the snake layout was measurably faster than the ladder layout. From this, we conclude that even a modest dependence on non-neighbor communication, such as in our ladder layouts, can have a detrimental effect on performance, and that snake layouts should be used when possible.

Platform Comparison

An additional observation in this study is the significant discrepancy in performance between the Raijin and Fujin platforms. Our estimate of cycle count indicates that most of this difference can be attributed to CPU clock speed. However, the x86 CPUs still demonstrate somewhat better performance per cycle. This can be attributed to various factors, such as differences in CPU cache architecture or compilers, but a more thorough investigation would be required to confirm the level of efficiency and underlying cause. An unexplored opportunity for improved hardware performance on Fujin is the Tofu barrier acceleration of the MPI reductions and broadcasts. Such operations are a dominant part of most diagnostic calculations, and the Tofu barrier could have an major impact on future simulations with high levels of diagnostic output.

Another point of consideration is the increasing trend of greater numbers of cores on each node, including MIC architectures such as Xeon Phi accelerators. Although Raijin outperformed Fujin in runtime, scalability was not possible without the introduction of redundant computational cores. No such additional resources were required to achieve scalable results on Fujin.

Submodel Performance

The difference in ocean scalability versus the ice and coupler submodel also indicates a potential target for optimization. Running the ocean and sea ice serially on each process is a potential bottleneck that could be addressed by moving the sea ice calculations onto separate core. In such a configuration, load balancing would become a greater concern, and a separation of ocean and ice would put a greater burden on the coupler. But given the inability of SIS to scale beyond 480 CPUs, and the very large number of sea ice timesteps required per ocean timestep, this may become a requirement in the future. The poor scalability of the coupling subroutines indicate that a parallelisation of of ocean and sea ice timesteps will not be sufficient to achieve future scalability. Efficient field exchange between models must also be a target of future optimisation efforts.

Conclusion

We have shown that the global 0.25° resolution ocean-sea ice configurations of the MOM ocean model can run efficiently on both x86 PRIMERGY and SPARC FX10 platforms. On the x86-based Raijin system, we were able to simulate over 10 model years per day at a sufficiently high level of efficiency. On the FX10 Fujin platform, we can simulate over 4 model years per day at a similar level of efficiency.

In order to achieve a high level of performance on the Raijin platform, we were required to either reduce the number of model processes per node to 12, leaving four cores to run idle during the simulation, or to enable hyperthreading on our CPUs, allowing two computational threads per core. Although reducing the number of cores yields a greater level of performance, hyperthreaded runs yield comparable results without any idle cores, thereby achieving much greater efficiency. When all 16 cores on each node were used without hyperthreading, the runtimes increased by a factor of three and were unable to scale beyond 480 CPUs.

Comparison of the submodels within the ice-ocean configuration shows that the ocean model scales up to 1920 cores with only a negligible loss of efficiency, while the sea ice and coupling components cease to scale after 480 CPUs. Future efforts to improve the scalability of MOM should target these subcomponents.

Both platforms have proven capable of running global high-resolution simulations of the ocean, and enable scientists to investigate the impacts of turbulent-scale processes on global climate. Development towards resolutions of 0.1° and beyond will require further investigations of the underlying codes and their resource requirements, and ongoing collaborations between industry and the academic research community.

Fig. 1 .

 1 Fig. 1. Node layouts for each Tofu unit. Each node number corresponds to 16 domain tiles.Figure (a) shows the snake layout and figure (b) shows the ladder layout.

  Fig. 1. Node layouts for each Tofu unit. Each node number corresponds to 16 domain tiles.Figure (a) shows the snake layout and figure (b) shows the ladder layout.

FigureFig. 2 .

 2 Figure 2b attempts to clarify this performance difference across the platforms by estimating the mean number of cycles per CPU with the following formula: 𝑁𝑁 𝑐𝑐 =f(1 -𝑝𝑝)𝜏𝜏 𝑜𝑜

Fig. 3 .

 3 Fig. 3. Model walltime (a) and scaling efficiency (b)

Fig. 4 .

 4 Fig. 4. Scaling efficiency for the ocean (a), sea ice (b), and coupler (c) submodels.

adfa, p. 1, 2011. © Springer-Verlag Berlin Heidelberg 2011

Acknowledgements

This work was conducted as part of the ACCESS Optimisation Project, a collaboration between NCI, Fujitsu, and the Australian Bureau of Meteorology. Research was undertaken with the assistance of resources from NCI, which is supported by the Australian Government. We are grateful to Mark Cheeseman of NCI, Symon Swann of Fujitsu Australia, and Motohiro Yamada of Fujitsu Japan for constructive reviews on an earlier draft of this paper.