WO2011094802A1 - Estimation of weather modification effects - Google Patents

Abstract

Disclosed is a method of estimating the effect of a weather modification apparatus over a trial area. The method comprises determining an observational model for meteorological observations at sites over the trial area using variables unrelated to the operation of the weather modification apparatus and variables related to the operation of the weather modification apparatus; determining a first set of meteorological values in the trial area using the observational model excluding the variables unrelated to the operation of the weather modification apparatus; and determining an estimate of the effect of the weather modification apparatus over the trial area using the first set of meteorological values.

Description

ESTIMATION OF WEATHER MODIFICATION EFFECTS

Technical Field of the Invention

The present invention relates generally to weather modification and, in particular, to measurement of the effects of weather modification technologies using statistical techniques.

Background

Weather modification technology is an area of increasing research and development activity, spurred in part by concern about climate change and its effect on local weather conditions. One example of particular interest to the agricultural industry is rainfall enhancement technology. However, investment in such technologies depends on conclusive demonstrations of their efficacy. In addition, assuming efficacy can be demonstrated, setting up a market for using and trading in such technologies depends on the ability to measure, or at least reliably estimate, their effects on a given area under given circumstances.

Weather systems are extremely complex, involving a great many parameters that can vary dramatically over vast volumes and time scales. Convincing theoretical explanations of the causal operation of current weather modification technologies have for this reason remained elusive. Computer simulations and laboratory models of the interaction of weather systems and weather modification technologies are inadequate as proving grounds for a variety of reasons. Demonstrating their efficacy and estimating their effects therefore rests, and will continue to rest, on statistical evaluation of data obtained while operating the technology under real-world conditions.

In previous trials of rainfall enhancement technologies, statisticians have relied on comparisons of trial results with long-term averages of rainfall on a given catchment. However, the high natural variability of rainfall data, both over time and over areas, and the influence of longer-term climate change has hindered conclusive demonstrations of efficacy using such techniques. There is a need for better statistical evaluation methods that can approach the conclusive evidence provided by real-time controlled laboratory experiments on other, smaller-scale technologies.

Summary

It is an object of the present invention to substantially overcome, or at least ameliorate, one or more disadvantages of existing arrangements.

Disclosed herein are methods for demonstrating whether, and estimating how much, rainfall enhancement technologies affect rainfall patterns in a given area. The disclosed methods use a correlation between observations of rainfall at different locations and / or at specified time intervals (24 hours is typical) to make concurrent predictions of rainfall in a target area. An advantage of the disclosed methods is that the methods do not simply attribute an observed change in the level of rainfall to an applied enhancement technology. Rather, the disclosed methods compare an explained variation in gauge to gauge rainfall with and without the enhancement technology in order to predict an increase in the probability and / or volume of rainfall due to operation of the enhancement technology. The disclosed methods are appropriate for typical datasets obtained during weather modification operations, due to the potential presence of spatio-temporal correlations in rainfall within and across the control and target areas. More generally, the disclosed methods allow more reliable and timely estimation of the effects of weather modification technologies, a result that is critical to the commercial application of such technologies. According to a first aspect of the present disclosure, there is provided a method of estimating the effect of a weather modification apparatus over a trial area, the method comprising:

determining an observational model for meteorological observations at sites over the trial area using variables unrelated to the operation of the weather modification apparatus and variables related to the operation of the weather modification apparatus; determining a first set of meteorological values in the trial area using the observational model excluding the variables unrelated to the operation of the weather modification apparatus; and

determining an estimate of the effect of the weather modification apparatus over the trial area using the first set of meteorological values.

According to a second aspect of the present disclosure, there is provided a method of estimating the effect of a weather modification apparatus over a trial area comprising at least one control area and a target area, the method comprising:

determining a control model by regressing meteorological observations in the trial area against variables unrelated to the operation of the weather modification apparatus; determining an effects model by regressing the meteorological observations in the trial area against the variables unrelated to the operation of the weather modification apparatus, and variables related to the operation of the weather modification apparatus; determining a first set of meteorological values in the target area using the control model and a second set of meteorological values in the target area using the effects model; and

determining a difference between the first set of meteorological values and the second set of meteorological values, the difference representing an estimate of the effect of the weather modification apparatus over the trial area. According to another aspect of the present disclosure, there is provided a computer program product including a computer readable medium having recorded thereon a computer program for implementing any one of the aforementioned methods.

Brief Description of the Drawings

At least one embodiment of the present invention will now be described with reference to the drawings, in which:

Fig. 1 is an illustration of an ground-based rainfall enhancement apparatus, according to one example;

Figs. 2A and 2B form a schematic block diagram of a general purpose computer system upon which the disclosed methods can be practised;

Fig. 3 shows a representation of an trial area comprising two control areas and a target area, according to one example;

Fig. 4 shows a map of a plume from the example apparatus illustrated in Fig. 1 ;

Fig. 5 shows an example dynamic partitioning of an area into target and control areas based on a steering wind direction;

Fig. 6 shows a further example dynamic partitioning of an area into target and control areas based on a steering wind direction;

Fig. 7 is a flow diagram illustrating a method of analysing rainfall and other meteorological observations obtained from rain gauges and stations within a trial area at specified analysis intervals over a trial period in order to estimate an increase in the level of rainfall within the trial area over the trial period due to the operation of a ground-based rainfall enhancement apparatus; Fig. 8 is a flow diagram illustrating a first method of determining a downwind direction;

Fig. 9 is a flow diagram illustrating a second method of determining a downwind direction;

Fig. 10 is a flow diagram illustrating an alternative method of analysing meteorological observation data obtained from sites within the trial area at the specified analysis intervals over the trial period;

Fig. 11 is a flow diagram illustrating a method of analysing rainfall and other meteorological observations obtained from rain gauges and stations within a trial area at specified analysis intervals over a trial period in order to estimate an increase in the probability that a rainfall event occurred during the trial period within the trial area due to the operation of a ground-based rainfall enhancement apparatus;

Fig. 12 is a flow diagram illustrating a method of determining a bias correction factor as used in the method of Fig. 10; and

Fig. 13 illustrates the radial tiling of the downwind region of a trial area into spatiotemporal groups of sites, as realised in a northwest steering wind.

Detailed Description including Best Mode

Where reference is made in any one or more of the accompanying drawings to steps and/or features, which have the same reference numerals, those steps and/or features have for the purposes of this description the same function(s) or operation(s), unless the contrary intention appears.

Methods 700, 1000, and 1 100 of analysing rainfall and other meteorological observations in order to estimate the effects of weather modification technology are described below with reference to Figs. 7, 10, and 1 1 respectively. The disclosed methods 700, 1000, and 1100 can be applied across a wide range of weather modification technologies, but are particularly suited to a ground-based rainfall enhancement apparatus (hereafter GREA) that emits an ion plume. The plume interacts with clouds downwind of the GREA to increase the probability, and subsequent intensity, of rainfall in a target area. This plume can include, but is not limited to, physical particulates such as silver iodide, or negative ions that attach to naturally occurring aerosols such as water vapour.

Fig. 1 is an illustration of a GREA 100 according to one example. The GREA 100 comprises a high voltage DC generator connected to a network of thin metal wire conductors supported by a framework 110 and surmounted by a set of pyramid-shaped structures 120. The typical approximate dimensions of the GREA 100 are 12 metres by 4 metres by 5 metres in height. The GREA 100 typically consumes about 500W of power and generates voltages of 80 to 85 kV.

The operation of the GREA 100 is as follows:

• The GREA 100 generates negative ions by corona discharge from the conductors into the surrounding air.

• The generated ions become attached, and transfer their electric charge, to atmospheric particles (aerosols).

• The charged aerosols form a plume which is conveyed to the upper atmosphere by wind, convection, and turbulence.

· The charged aerosols influence the collision and coalescence of cloud droplets in the upper atmosphere, depending on the magnitude of the charge.

• The coalesced droplets precipitate as rain downwind from the GREA 100.

Figs. 2A and 2B collectively form a schematic block diagram of a general purpose computer system 200, upon which the methods 700, 1000, and 1100 can be practised. As seen in Fig. 2A, the computer system 200 is formed by a computer module 201, input devices such as a keyboard 202, a mouse pointer device 203, a scanner 226, a camera 227, and a microphone 280, and output devices including a printer 215, a display device 214 and loudspeakers 217. An external Modulator- Demodulator (Modem) transceiver device 216 may be used by the computer module 201 for communicating to and from a communications network 220 via a connection 221. The network 220 may be a wide-area network (WAN), such as the Internet or a private WAN. Where the connection 221 is a telephone line, the modem 216 may be a traditional "dial- up" modem. Alternatively, where the connection 221 is a high capacity (eg: cable) connection, the modem 216 may be a broadband modem. A wireless modem may also be used for wireless connection to the network 220.

The computer module 201 typically includes at least one processor unit 205, and a memory unit 206 for example formed from semiconductor random access memory (RAM) and semiconductor read only memory (ROM). The module 201 also includes an number of input/output (I/O) interfaces including an audio-video interface 207 that couples to the video display 214, loudspeakers 217 and microphone 280, an I/O interface 213 for the keyboard 202, mouse 203, scanner 226, camera 227 and optionally a joystick (not illustrated), and an interface 208 for the external modem 216 and printer 215. In some implementations, the modem 216 may be incorporated within the computer module 201, for example within the interface 208. The computer module 201 also has a local network interface 21 1 which, via a connection 223, permits coupling of the computer system 200 to a local computer network 222, known as a Local Area Network (LAN). As also illustrated, the local network 222 may also couple to the wide network 220 via a connection 224, which would typically include a so-called "firewall" device or device of similar functionality. The interface 21 1 may be formed by an Ethernet™ circuit card, a Bluetooth™ wireless arrangement or an IEEE 802.1 1 wireless arrangement.

The interfaces 208 and 213 may afford either or both of serial and parallel connectivity, the former typically being implemented according to the Universal Serial Bus (USB) standards and having corresponding USB connectors (not illustrated). Storage devices 209 are provided and typically include a hard disk drive (HDD) 210. Other storage devices such as a floppy disk drive and a magnetic tape drive (not illustrated) may also be used. An optical disk drive 212 is typically provided to act as a non- volatile source of data. Portable memory devices, such optical disks (eg: CD-ROM, DVD), USB-RAM, and floppy disks for example may then be used as appropriate sources of data to the system 200.

The components 205 to 213 of the computer module 201 typically communicate via an interconnected bus 204 and in a manner which results in a conventional mode of operation of the computer system 200 known to those in the relevant art. Examples of computers on which the described arrangements can be practised include IBM-PC's and compatibles, Sun Sparcstations, Apple Mac™ or alike computer systems evolved therefrom.

The methods 700, 1000, and 1100 may be implemented using the computer system 200 as one or more software application programs 233 executable within the computer system 200. In particular, the steps of the methods 700, 1000, and 1 100 are effected by instructions 231 in the software 233 that are carried out within the computer system 200. The software instructions 231 may be formed as one or more code modules, each for performing one or more particular tasks. The software may also be divided into two separate parts, in which a first part and the corresponding code modules performs the methods 700, 1000, and 1100 and a second part and the corresponding code modules manage a user interface between the first part and the user.

The software 233 is generally loaded into the computer system 200 from a computer readable medium, and is then typically stored in the HDD 210, as illustrated in Fig. 2 A, or the memory 206, after which the software 233 can be executed by the computer system 200. In some instances, the application programs 233 may be supplied to the user encoded on one or more CD-ROM 225 and read via the corresponding drive 212 prior to storage in the memory 210 or 206. Alternatively the software 233 may be read by the computer system 200 from the networks 220 or 222 or loaded into the computer system 200 from other computer readable media. Computer readable storage media refers to any storage medium that participates in providing instructions and/or data to the computer system 200 for execution. and/or processing. Examples of such storage media include floppy disks, magnetic tape, CD-ROM, a hard disk drive, a ROM or integrated circuit, USB memory, a magneto-optical disk, or a computer readable card such as a PCMCIA card and the like, whether or not such devices are internal or external of the computer module 201. Examples of computer readable transmission media that may also participate in the provision of software, application programs, instructions and/or data to the computer module 201 include radio or infra-red transmission channels as well as a network connection to another computer or networked device, and the Internet or Intranets including e-mail transmissions and information recorded on Websites and the like.

The second part of the application programs 233 and the corresponding code modules mentioned above may be executed to implement one or more graphical user interfaces (GUIs) to be rendered or otherwise represented upon the display 214. Through manipulation of typically the keyboard 202 and the mouse 203, a user of the computer system 200 and the application may manipulate the interface in a functionally adaptable manner to provide controlling commands and/or input to the applications associated with the GUI(s). Other forms of functionally adaptable user interfaces may also be implemented, such as an audio interface utilizing speech prompts output via the loudspeakers 217 and user voice commands input via the microphone 280.

Fig. 2B is a detailed schematic block diagram of the processor 205 and a

"memory" 234. The memory 234 represents a logical aggregation of all the memory devices (including the HDD 210 and semiconductor memory 206) that can be accessed by the computer module 201 in Fig. 2A.

When the computer module 201 is initially powered up, a power-on self-test (POST) program 250 executes. The POST program 250 is typically stored in a ROM 249 of the semiconductor memory 206. A program permanently stored in a hardware device such as the ROM 249 is sometimes referred to as firmware. The POST program 250 examines hardware within the computer module 201 to ensure proper functioning, and typically checks the processor 205, the memory (209, 206), and a basic input-output systems software (BIOS) module 251, also typically stored in the ROM 249, for correct operation. Once the POST program 250 has run successfully, the BIOS 251 activates the hard disk drive 210. Activation of the hard disk drive 210 causes a bootstrap loader program 252 that is resident on the hard disk drive 210 to execute via the processor 205. This loads an operating system 253 into the RAM memory 206 upon which the operating system 253 commences operation. The operating system 253 is a system level application, executable by the processor 205, to fulfil various high level functions, including processor management, memory management, device management, storage management, software application interface, and generic user interface.

The operating system 253 manages the memory (209, 206) in order to ensure that each process or application running on the computer module 201 has sufficient memory in which to execute without colliding with memory allocated to another process. Furthermore, the different types of memory available in the system 200 must be used properly so that each process can run effectively. Accordingly, the aggregated memory 234 is not intended to illustrate how particular segments of memory are allocated (unless otherwise stated), but rather to provide a general view of the memory accessible by the computer system 200 and how such is used.

The processor 205 includes a number of functional modules including a control unit 239, an arithmetic logic unit (ALU) 240, and a local or internal memory 248, sometimes called a cache memory. The cache memory 248 typically includes a number of storage registers 244 - 246 in a register section. One or more internal buses 241 functionally interconnect these functional modules. The processor 205 typically also has one or more interfaces 242 for communicating with external devices via the_ system bus 204, using a connection 218.

The application program 233 includes a sequence of instructions 231 that may include conditional branch and loop instructions. The program 233 may also include data 232 which is used in execution of the program 233. The instructions 231 and the data 232 are stored in memory locations 228-230 and 235-237 respectively. Depending upon the relative size of the instructions 231 and the memory locations 228-230, a particular instruction may be stored in a single memory location as depicted by the instruction shown in the memory location 230. Alternately, an instruction may be segmented into a number of parts each of which is stored in a separate memory location, as depicted by the instruction segments shown in the memory locations 228-229.

In general, the processor 205 is given a set of instructions which are executed therein. The processor 205 then waits for a subsequent input, to which it reacts to by executing another set of instructions. Each input may be provided from one or more of a number of sources, including data generated by one or more of the input devices 202, 203, data received from an external source across one of the networks 220, 222, data retrieved from one of the storage devices 206, 209 or data retrieved from a storage medium 225 inserted into the corresponding reader 212. The execution of a set of the instructions may in some cases result in output of data. Execution may also involve storing data or variables to the memory 234.

The disclosed methods use input variables 254, that are stored in the memory 234 in corresponding memory locations 255-258. The disclosed methods produce output variables 261 , that are stored in the memory 234 in corresponding memory locations 262- 265. Intermediate variables may be stored in memory locations 259, 260, 266 and 267.

The register section 244-246, the arithmetic logic unit (ALU) 240, and the control unit 239 of the processor 205 work together to perform sequences of micro-operations needed to perform "fetch, decode, and execute" cycles for every instruction in the instruction set making up the program 233. Each fetch, decode, and execute cycle comprises:

(a) a fetch operation, which fetches or reads an instruction 231 from a memory location 228;

(b) a decode operation in which the control unit 239 determines which instruction has been fetched; and

(c) an execute operation in which the control unit 239 and/or the ALU 240 execute the instruction.

Thereafter, a further fetch, decode, and execute cycle for the next instruction may be executed. Similarly, a store cycle may be performed by which the control unit 239 stores or writes a value to a memory location 232. Each step in the methods 700, 1000, and 1100 is associated with one or more segments of the program 233, and is performed by the register section 244-247, the ALU 240, and the control unit 239 in the processor 205 working together to perform the fetch, decode, and execute cycles for every instruction in the instruction set for the noted segments of the program 233.

The methods 700, 1000, and 1100 may alternatively be implemented in dedicated hardware such as one or more integrated circuits performing the functions or sub functions of Fig. 7. Such dedicated hardware may include graphic processors, digital signal processors, or one or more microprocessors and associated memories.

To collect data on which to carry out the analysis methods 700, 1000, and 1100 described below, firstly a trial area is chosen. The trial area can range from hundreds to thousands of square kilometres. Typically, a trial- area is chosen because it would benefit from an increase in rainfall, and is under the influence of a single atmospheric air mass.

The methods 700, 1000, and 1 100 use data, in the form of historical records, from a plurality (typically several dozen) of meteorological stations and rain gauges (collectively referred to as sites) within the designated trial area. The amount of sites preferred typically depends on the spatial variation in meteorology and in particular rainfall, across the trial area. In a trial area of steep and undulating topography, it is preferable to have a higher number of sites at higher spatial density to account for topographic induced variations in meteorology and rainfall over small scales. While in a meteorologically homogenous trial area, a smaller number of sites relative to the area may be sufficient. In any case however, the network of sites needs to be large enough to adequately measure natural rainfall variations within the trial area. A minimum of thirty (30) years of historical data is preferable, however, periods shorter than this may be used if too few long-term stations are to be found within the trial area. Typically, the most comprehensive weather records are maintained by government meteorological services. However, data can be obtained from private individuals or companies conducting their own meteorological observations if it is possible to conduct quality control and verify the historical data.

When assessing the suitability of a trial area, the number and distribution of sites is considered. The more sites and the denser the network thereof, the better suited that trial area will be to the analysis methods 700, 1000, and 1100 described below. In instances where existing sites are sparse, additional rain gauges may be placed in the trial area. The placement of gauges has an effect on the analysis, and a statistical technique is applied to the historical rainfall and wind data from the existing sites to account for the spatial dispersion pattern of the gauges. Typically, rain gauges are set to record rainfall at a minimum of twenty-four (24) hourly intervals (usually from 9 a.m. to 9 a.m.). However, optimal data recording is at ten (10) minute intervals, which then can be averaged or aggregated to get hourly or daily data.

Once a trial area has been chosen, target and control areas may be assigned within the trial area. The target area is nominally under the influence of the GREA, while the control area(s) are nominally free of the influence of the GREA. The control and target areas are also preferably highly correlated in rainfall at a specified analysis interval. To confirm this, correlation analysis is conducted using data from sites within each proposed control and target area. Spatial correlation tends to decline with distance. The closer together a well-defined target and control area are, the greater the level of experimental control. Exploiting the directional aspect of a GREA within the target area allows the target area to be partitioned into segments or sub-areas where the enhancement effect should be greater or less, within relatively close proximity to the GREA. Fig. 3 shows a representation of an example trial area 300 comprising two control areas (north 320 and south 330), and a target area 310 surrounding a GREA 305, of which an example is the GREA 100 of Fig. 1. A high density network of sites, e.g. 340, is shown within each area.

Spatial correlation can vary over different time intervals, giving rise to spatio- temporal correlation. For example, at an hourly time interval, the spatial correlation between observed rainfalls may be much lower than the spatial correlation that might be observed over a monthly or annual time interval. Again, for example spatial correlation may be affected by wind direction and speed and the latter of the two changing over time. The amount of non-zero spatial correlation between rainfall values at the chosen analysis interval will determine how accurately measured rainfall can be used to estimate the change in rainfall attributed to the GREA 305. That is, spatio-temporal correlation determines in large part the level of real time experimental control that can be achieved in a rainfall enhancement experiment. However, the effectiveness of this control depends on the level of spatial correlation. A typical analysis interval at which spatial correlation is sufficient to provide valuable information is one day (24 hours).

Using the historical rainfall and wind data from the stations 340, an analysis is conducted to determine the "best" position for the GREA 305 to affect the chosen target area. The best position for the GREA 305 is considered to be a position that is upwind of the target area when the majority of precipitation occurs. For example, if the majority of precipitation occurs when the wind is from the southwest, then the best position for the GREA 305 would be to the southwest of the target area. In the example of Fig. 3, there is no predominant wind direction during periods of precipitation, so the best position for the GREA 305 is in the centre of the target area 310.

Variations in topography, area of need, client requirements, or variability in weather patterns may necessitate the use of more than one GREA 305. The analysis methods 700, 1000, and 1 100 described below can be applied to trials using more than one GREA 305 so long as the downwind area of effect for each individual GREA 305 can be defined. In some cases it may be necessary to aggregate these areas, when the individual areas of effect overlap. However, for the purpose of this disclosure one GREA 305 will be assumed.

Typically, air masses can be parametrised by their wind speed and direction, humidity, and pressure at various levels in the atmosphere. Wind speed and direction are measured at a number of stations (typically greater than 3, so as to obtain an average local correlation) closest to the location of the GREA 305 and at corresponding relative locations in the control area(s). Typically the chosen stations are within a 10km range of the GREA 305, but this range will be dependent on local station distribution. To conduct the wind direction analysis as described below, the measured wind directions are typically binned into one of eight compass bearings, those being N, NE, E, SE, S, SW, W and NW. Suitable correlation coefficients vary from location to location.

Rainfall observation records are taken at a number (typically greater than 3, so as to obtain an average local correlation) of sites (e.g. 340) closest to the location of the GREA 305 and at corresponding relative locations in the control area(s) 320 and 330. Typically the chosen stations are within a 10 kilometre range of the GREA 305, but this range will be dependent on local station distribution. To conduct the rainfall analysis as described below, the rainfall observation records are typically taken at monthly and daily intervals during the trial period.

The potential rainfall enhancement effect of a GREA 305 may be separated into two components, as follows:

(1) A change in the probability that within a specified interval, typically 24 hours, a rainfall event will occur; (2) Given that a rainfall event does occur within that interval, a change in the expected level of rainfall over that interval.

The methods 700 and 1000 described below estimate only the second of these two components. The method 1 100, described below with reference to Fig. 1 1, estimates the first of these two components.

The methods 700, 1000, and 1 100 described below take into account the spatio- temporal correlation between rainfall observations at different locations at the specified analysis interval. The methods 700, 1000, and 1100 also take into account the fact that the direction and location of any enhancement effect of a GREA (e.g. the GREA 100) depends on:

(1) Meteorological conditions, such as wind direction and speed;

(2) Fixed effects of topography and land cover giving rise to orographic - lifting and affecting turbulence.

A preparatory step for the methods 700, 1000, and 1 100 described below is to determine the downwind direction in the vicinity of the GREA 305, in order to determine the area of effect, which in essence is a dynamic target area, continuously defined by winds within the broader target area 310.

For the determination of downwind direction, it is necessary to observe the direction and speed of surface and upper level winds, as these can vary with height. To eliminate any possible bias, either of two methods 800 (see Fig. 8) and 900 (see Fig. 9) is used to determine the downwind direction and hence define the dynamic target area. The first method 800 uses wind observation data from local stations and radiosonde data, typically from weather balloons. The downwind direction so determined is referred to herein as the steering wind direction. The second method 900 additionally uses numerical modelling. The downwind direction determined by the second method 900 is referred to herein as the principal direction.

One or more steps of the method 800 may be implemented as software resident on the hard disk drive 210 and being controlled in its execution by the processor 205. In the first step 805 of the method 800, the processor 205 obtains weather observations from the stations (e.g. 340) in the target area 310, including wind speed and direction, temperature, pressure and humidity. Data values representing these weather observations may be obtained and stored in the memory 206 of the computer system 200. In one implementation the weather observation data values may be downloaded by the processor 205 over the network 220, where each of the stations (e.g. 340) is linked to the network 220. At the next step 810, the wind observation data is processed by the processor 205 along with other trial area wind observations to determine a spatially and temporally averaged surface wind speed and direction for the specified analysis interval, typically one day. At the next step 815 in the method 800, the processor 205 determines a representative vertical wind profile using radiosonde observation data. A typical wind profile comprises wind speeds and directions measured at various heights in the atmosphere. Common radiosonde data includes wind speed and direction at standard pressure levels, as well as other levels. A pressure level is a level of defined pressure (ex: 300 hPa) and the height (above MSL) at which that pressure value is found. It is typically used because it is easier for balloon sounding using pressure sensors to note pressure(s) rather than the height(s) at which the data is recorded. The standard pressure levels typically include 850 hPa, 700 hPa, 500 hPa. However modern measuring equipment (such as acoustic sounders) measure wind speed and direction directly at heights and can be readily used.

Two options exist to obtain this data. If an existing radiosonde balloon is launched in the vicinity of the GREA 305, and is considered representative of the trial area air mass, then its data can be used. Alternatively, it may be necessary to launch a radiosonde balloon from the trial area to obtain the observation data.

Once the representative wind profile is obtained and stored in the memory 206, the steering wind direction is determined by the processor 205 at step 820 as a speed weighted average of the representative wind profile in a layer of the atmosphere considered to contain the majority of the dispersion of the ion plume, including the layers where the cloud occurs. This layer varies from location to location, but typically in tropical regions is the lower 6 kilometres (20,000 feet) of the atmosphere, or below the 500 hPa level, while in more mid-latitude areas the layer may be taken as below (10,000 feet) or 700 hPa level. In areas of complex topography where surface observations may be highly variable, the layer may be defined with a lower bound as well, for example between the 850 hPa to 700 hPa levels.

One or more steps of the method 900 of determining the principal direction may also be implemented as software resident on the hard disk drive 210 and being controlled in its execution by the processor 205. In the first step 905 of the method 900, the processor 205 obtains weather observations from the stations (e.g. 340) in the target area 310, including wind speed and direction, temperature, pressure and humidity, as in step 805 of the method 800. The processor 205 then uses the weather observation data as input for numerical meteorology models to determine the principal direction of the plume generated by the GREA 305. A number of atmospheric dispersion models could be used for this task, including, but not limited to, ADMS 3, ATSTEP, AUSTAL2000, CALPUFF, DISPERSION21, ISC3, MEMO Model, MERCURE, PUFF-PLUME, RIMPUFF, and SAFE_AIR. These models can be used as a mathematical simulation (of varying complexity) of how particulates, which are used as a proxy for ions generated from the GREA, disperse in the ambient atmosphere. The simulation is performed with computer programs that solve the mathematical equations and algorithms which simulate the particle dispersion. The dispersion models are used to estimate or to predict the downwind concentration and thus location, of the ions (and resulting charged aerosols) generated from the GREA. However the majority of dispersion models are designed to predict pollutant dispersion and ground based pollutant concentrations so do not include simulation of atmospheric processes above the boundary layer such as deep convection. Preferably, full three-dimensional mesoscale meteorology modelling suites such as MM5 or WRF are run independently to obtain the most accurate and relevant results for the local area.

Once the chosen atmospheric dispersion or full meteorology model output has completed by the processor 205 in step 910 of the method 900, yielding predictions of ion concentrations (and ion plume location) emitted frorn the GREA at each time step, the processor 205 in step 915 determines the principal direction of the modelled plume. Using the 1 -dimensional and simple 2D models, the principal plume direction is an output of the model. Alternatively in more complex 3D full meteorology modelling suites, the processor 205 determines the principal direction from the ion concentration (plume) data. In one implementation, this can be depicted graphically and is typically done by drawing two lines outward from the GREA and which tangent the ion plume and extend to the limit of the plume, so creating a bounding triangle. The principal plume direction is then calculated as the direction, within this triangle, such that equal area of ion plume is located either side of the line drawn from the GREA in that direction.

Fig. 4 shows the map of a MM5-modelled plume 400 from the GREA 100 illustrated in Fig. 1. The tangent lines 410 and 420 and resultant principal plume direction are also shown. For this output the lower level winds were from the WNW and the mid to upper level winds were from the NW. The resultant principal direction is shown as the middle dashed line 430. The production of Fig 4, a MM5 model run, used a modelling domain of 0.5S x 0.65E centred on 35.27S, 138.75S, with a spatial resolution of 0.005 deg in the x- and y-directions and 12 levels in the vertical. The model results were output as 10 minute averages. A number of assumptions were made in the modelling study: (1) Ions were treated as passive tracers. This means that: (a) ions were assumed to behave the same as gaseous pollutants or particles; (b) recombination of ions in the atmosphere was neglected; (c) deposition of ions (including deposition through rainfall) was also neglected. Full reflection is assumed for particles that intersect the ground or model-top. The GREA was modelled as area sources with a height of 5 metres above ground level and a surface area of 36 square metres.

The use of full 3D meteorology modelling suites yields more accurate and detailed depiction of the ion plume, particularly in areas of complex terrain and meteorology. However these suites are commonly expensive and time consuming to run, and so the use of steering wind direction as the downwind direction is preferred when the output is considered an adequately accurate representation of the advection of the ion plume in the atmosphere, a consideration which will vary between each location and time of year.

The target area 310 is partitioned into wind flow sectors using the downwind direction determined using the first method 800 or the second method 900. Partitioning the target area 310 serves to focus the signal generated by the GREA 305 and allows spatio- temporal correlation structure to be identified and measured. The partitioning is effectively dynamic because the downwind direction varies between analysis intervals. As a consequence, the fact that the GREA 305 may be located near land features that enhance or limit rainfall becomes a small source of potential bias. Sectors are defined by a set of angles, increasing symmetrically in absolute terms about the downwind direction. Each site is then classified according to the sector containing its location. In one example as shown in Fig. 5, the circular target area 310 is partitioned into five (5) sectors, based on an average downwind direction that is due northerly (0 degrees) over the trial period.

Sector 1 : 0° < ± angle < 30 °

Sector 2: 30° < ± angle < 60 °

Sector 3 : 60° < ± angle < 90 °

Sector 4: 90° < ± angle < 135 °

Sector 5: 135° < ± angle < 180 °

As illustrated in Fig. 6, Sectors 1 and 2 form a dynamic target area 600 labelled as downwind, sectors 3 and 4 are labelled as crosswind, and sector 5 is labelled as upwind.

In a second implementation, a site is considered to be downwind for an analysis interval when both the following conditions are satisfied: sin(<9 - co)(lat - lat_A ) + cos(# - <y)(/on - long_A ) < 0;

sin(# + co)(lat - lat_A ) - cos(d? + o)[long - long_A ) < 0. where Θ is the angle from the wind direction vector (in the case of the downwind sector, this value is 30°), ω is the bearing of the site from the GREA location (0^° being due north), and lat and long and lat_A and long_A are the latitude and longitude of the site and the

GREA respectively. Any sites not within the downwind sector are considered to be crosswind/upwind.

When there are two GREA sites in use during a trial, site observations during an analysis interval may be classified according to their downwind and crosswind / upwind orientations relative to the two GREA locations. The definition of these "spatio-temporal" classifications is aimed at grouping downwind sites that are likely to have had similar exposure to prevailing meteorological conditions, including exposure to a GREA. The classifications are defined based on the radial angle (C26) made by a site with the first GREA (C2) at a first location relative to the average steering wind direction for the analysis interval, and the same angle (C3<9) relative to the second GREA C3) at a second location. Six classifications are defined as listed in Table 1 below:

; _C20 greater than 30° and C36> greater than 60° Γ ^{' "}Ί ΐ^^~ ^ΞΞ] b C20 less than or equal to 30° and C30 greater than 60°

ί c C20 less than C36> and both C20, C30 less than or equal to ^'.^^≤^ ^' Γ Ι ^ΞΖΖΙ,^ " d C30 less than or equal to C20 and both C20, C30 less than or equal to 60°

: · e C3# less than or equal to 30°^'and C20 greater than 60°

f C30 greater than 30° and C20 greater than 60°

Table 1: Spatio-temporal classifications

The grouping of sites in accordance with the definitions in Table 1 effectively forms a radial tiling of the downwind region for each analysis interval, as illustrated in Fig.

13, which illustrates the radial tiling 1300 of the downwind region of a trial area into spatiotemporal classifications a to f as defined in Table 1 and realised in a northwest steering wind.

Rainfall sites in the target and control areas may also be dynamically classified based on other meteorological conditions besides wind direction. The dynamic site classification scheme addresses two key problems in estimating the amount of rainfall enhancement. The first is that climatic patterns over small geographic areas, such as a trial area, are not stable over time. It is difficult to make meaningful comparisons of how much it rained during the trial or trials to corresponding periods in the past, as there is a very high level of natural variation over time. Second, rainfall patterns are not geographically stable. There is a high level of natural variability in rainfall within geographic regions that tends to increase with the distance between regions. This again makes meaningful comparison between the trial area and other areas difficult. The dynamic classification scheme allows the use of nearby sites to obtain contemporaneous control and effect measurements. With multiple GREA sites and randomised cross-over experimental designs this may be limited to the area over which the downwind areas do not overlap. Given substantial areas without cross-overlap, the use of multiple sites that are near to each other as well as to the GREA site affords a better 'like with like' comparisons due to spatial correlation in daily rainfall patterns within a small geographic area.

The use of a dynamic classification also tends to result in having individual sites that provide control and effect measurements at different points in time over the trial. This reduces the likelihood that rainfall measurements will be biased due to factors related to site location, as for example being to the east or the west of the GREA location.

The classification scheme itself is based upon an evaluation of a range of factors including the distribution of wind directions throughout the atmosphere and the mathematical modelling of a point source plume of inert particles, as is often used to model point source pollution into the atmosphere.

Having winds from different directions over time during the trial reduces the likelihood that a given wind flow sector will be associated with any fixed orographic effects. For example, if rainfall was consistently associated with one wind direction then, for example, mountain ranges or hills could produce areas of elevated or reduced rainfall that were linked with a particular wind flow sector. It is preferable to have a wide distribution of downwind directions, so there is little scope for orographic bias. Moreover, having winds from different directions over time reduces the likelihood that a given wind flow sector will be associated with any fixed orographic effects.

Fig. 7 is a flow diagram illustrating the method 700 of analysing meteorological observation data obtained from the rain gauges and stations (e.g. 340) within the trial area (e.g. 300) at the specified analysis intervals over the trial period in order to estimate an increase in the level of rainfall within the trial area over the trial period due to the operation of a ground-based rainfall enhancement apparatus (e.g. 305). The method 700 may be implemented as software resident on the hard disk drive 210 and being controlled in its execution by the processor 205 of the computer system 200.

As described below, the analysis method 700 is specific to the purpose of determining whether the operation of rainfall enhancement technology (e.g. the GREA 100) is associated with increased rainfall in the target area (e.g. 310) or the dynamic target area (e.g. 600), conditional on the fact that a rainfall event has occurred, and if so, estimating the amount of increase. A rainfall event is defined as having at least one nonzero recorded measurement of rainfall within the trial area within the specified analysis interval. However, the method 700 is readily generalisable to estimating the weather . modification effects of other weather modification technologies.

The method 700 is described below with reference to a first trial of the GREA 100. Some of the main data elements of the first trial are summarised below:

Trial period: five months (144 days)

Analysis interval: 24 hours (9 a.m. to 9 a.m.)

A target area 310 was identified and defined as a 70km radius circle, centred on the GREA south west of Bundaberg, Queensland, (at 25° 21 ' 38.72"S, 151 ° 55' 15.19Έ) Two control areas: (1) northern control area 320 near Gladstone and (2) southern control area 330 near Gympie

Rainfall measurements from a total of 165 sites for the trial period:

Target Area: 1 17 gauges, South Control Area: 24 gauges, North Control Area: 24

Gauges

Rainfall events: 17499 site-interval records

Of these there were 5354 rainfall events (rainfall > 0 at a site in target or control areas). The distribution of events as follows: Target 3823, South control 831 , North Control 700.

Historical rainfall data available 1978 - 2007

Sites:

Elevation

Latitude and longitude

In the first step 705 of the method 700, the processor 205 determines a trend in observation data associated with rainfall events in the control areas. The trend is determined using a penalized spline that is fitted to the natural logarithm of the rainfall observations (LogRain) associated with rainfall events at each site, in order to reduce the relative influence of extreme rainfall events. The penalty factor is chosen in order to provide a compromise between a straight-line trend (high penalty factor) that ignores day- to-day variation in rainfall at different sites and an excessively variable trend (low penalty factor) that is completely driven by rainfall observation at sites in the control areas. The chosen penalty factor also removes some of the variability induced by intervals where only a few isolated non-zero rainfall observations are recorded.

In the second step 710, the processor 205 calculates the deviation in observed rainfall at each site in the control and target areas from the trend and regresses the resulting deviations against the latitude and longitude of the sites, an interaction of latitude and longitude effect, and a fixed interval or rainfall event effect, yielding a first regression model. The inclusion of the fixed interval or rainfall event effect accounts for short-term temporal variations about the trend. The fit of the first regression model gives an estimate of its explanatory power. In the first trial, the fit of the first regression model was 27.7%. The processor 205 then, at step 715, determines whether the explanatory power of the first regression model is satisfactory. It is not possible to measure all the characteristics that may influence rainfall or to measure them satisfactorily in every case. In particular, upper air measurements of temperature, humidity, pressure etc. are not reflected in the model except in so far as they may be correlated with the variables which are measured. However, the characteristics used in this first regression model generally explain around 28 per cent of the variation in rainfall that is observed. This indicates that 28% of rainfall variation is explained by the independent variables in the model. The remaining 72% consists of factors that were not measured. In research of highly variable systems, particularly where there are a large number of variables not explored in the research, a model which explains 40% of the variability or more is considered strong. For the particular studies conducted here, explanatory power of the order of 25% is generally held to be useful. In the first trial, the explanatory power, as estimated using the model fit, was found to be satisfactory.

If not, the processor 205 returns to step 710 to again determine the first regression model using a different combination of non-GREA related covariates. This iterative refinement of the first regression model continues until the explanatory power of the first regression model is satisfactory, when the method 700 proceeds to step 720.

In the following step 720, the processor 205 calculates expected rainfall at each site in the control and target areas at the analysis interval as the sum of a trend value at the site obtained from the spline fitted in step 705 and an estimate of the site deviation obtained from the first regression model computed at step 710. The differences between the control and target area deviations are attributable to location.

In an optional step 730, shown in dashed outline in Fig. 7, the processor 205 determines a second regression model in similar fashion to step 710, except including GREA-related variables, which typically include variables describing whether the GREA was operating and the amount of time it was operating during the analysis interval or previous analysis intervals (known as lags), the distance each site is from the GREA, and the wind flow sectors, as well as their interactions, as described in more detail below with reference to step 760. If the GREA-related variables do not influence rainfall, one would not expect the explanatory power of the non-GREA -related variables in the first regression model to be reduced, or diluted, by the inclusion of the GREA-related variables. Conversely, if the second regression model, including GREA-related variables, determined at step 730 has greater explanatory power than the first regression model, which includes only non-GREA -related variables, this is evidence of GREA influence on rainfall. Step 730 is performed for two approaches to defining which sites are controls, and which are GREA-influenced (i.e. in the target area), i.e. static (independent of wind direction) and dynamic (i.e. based on downwind direction as described above).

In a second optional step 735, shown in dashed outline in Fig. 7, the processor 205 compares the explanatory power of the second regression model, estimated as the model fit, with that of the first regression model. If the explanatory power of the second regression model is not greater, the method 700 halts at step 740. Otherwise, with some confidence that the GREA 100 does indeed have an effect on observed rainfall, the method 700 continues at step 750. In the first trial, the fit of the second regression model was found to be greater than the fit of the first regression model.

At step 750, the processor 205 determines a relationship known as the Control Model between non-zero rainfall observations at the individual sites in the control and target areas without taking into account any GREA effects, by regressing those observations against the expected values for the respective sites computed at step 720. The covariates in the Control Model also include historical average rainfall data, preferably comprising at least 30 years of monthly average rainfall data for the control areas as well as the target areas. However, shorter periods can be used. In absolute terms, these averages may not be reflective of the seasonal conditions that prevailed in the study regions. However, their relative levels may still explain regional differences in rainfall due to fixed effects such as topography. In the first trial, the fit of the Control Model was 38.5%.

At the next step 760, the processor 205 determines, in similar fashion to step 750, a model known as the Effects Model, which includes all the variables used in the Control Model, plus variables related to the operation of the GREA 305. In the first trial, the 'GREA-related variables' included in the Effects Model were:

• GREA Active (equal to one if the GREA was operating during the current analysis interval, or the two previous intervals)

• GREA Distance (Euclidean distance from the GREA site)

• Duration (minutes) of GREA operation during the analysis interval as well as for the previous analysis intervals (GREA Lags 0 to 2)

• GREA Distance * Duration interactions

• Wind flow Sector: Indicates site's angular distance from being directly downwind of the GREA during the analysis interval. The sectors are numbered from one to eight as described above with reference to Fig. 5.

• Lagged Wind flow Sector (previous analysis interval's values of Wind flow Sector)

• Wind flow Sector * Lagged Wind flow Sector interactions

Other GREA-related variables may be contemplated in step 760. In the first trial, using the above variables, the fit of the Effects Model was 40.0%.

In the next step 770, the processor 205 estimates the effect of the GREA by determining two values of (i.e. predicting) rainfall at each site in the target area at each analysis interval: one using the Effects Model determined in step 760 and the other using the Control Model determined in step 750. Since the two models model LogRain, the difference between the two predictions is also a logarithm. On back-transformation, with a nonparametric correction for the resulting transformation bias, the difference becomes a ratio representing the relative change in expected rainfall with and without the operation of the GREA. The ratio therefore represents the estimated enhancement effect of the GREA at each site and analysis interval. This approach has the advantage over conventional comparisons between predicted and observed rainfall, in that the difference between the models can only be attributed to the technology-related variables, i.e. the operation of the GREA itself. In conventional analysis methods, variation between the model predictions and the rainfall observations is likely to be wrongly attributed to the operation of the GREA.

Step 780 follows, at which the processor 205 estimates the on-ground environmental effect of the GREA. Using the enhancement effect at each site-interval estimated at step 770, an estimate of the volume of rainfall attributable to the GREA may be determined as the contribution to the total observed rainfall at a site that was due to the operation of the GREA:

Estimated GREA Contribution = (Estimated Enhancement Effect - 1) * Rainfall

The GREA contributions at a site may be aggregated over the full trial period or a part thereof.

The enhancement values estimated at step 770 can also be used to estimate the in- reservoir effect of the GREA. This is done for a given reservoir by applying local run-off coefficients to relevant site observations to determine the amount of rainfall attributable to the GREA which has entered the reservoir. A correction factor accounting for rainfall which fell directly into the reservoir needs to be applied, as this amount has a 100 per cent run-off coefficient (ignoring evaporation).

At the final step 790 of the method 700, the processor 205 calculates a standard error on the estimated percentage increase in rainfall due to the operation of the GREA. The variable %EI is an estimate of the percentage of rainfall in the sites in the area affected by the GREA that is due to the operation of the GREA. The processor 205 first determines %EI as a rainfall-weighted average of the enhancement effects estimated at step 770 for each record i (typically site-analysis interval), expressed as a percentage:

∑(£, - !) ,

%EI = 100 x -i-^ = Ι ΟΟ χ ΕΡΙ

i where ?, is the observed rainfall, E, is the Estimated Enhancement Effect, EPI is the estimated proportional effect, and the summation is over all records of interest, e.g. all site observation data over the trial period, or site observation data for a particular month and for a particular Wind flow Sector, or for their cross-classification (classification according to both Wind Flow Sector and Month at the same time). The processor 205 determines the standard deviation of %EI via Taylor series linearization as follows: σ_ε (%EI) = 100 x a_e (EPI) = 100 x i x a_e (z )

R where R is the average rainfall for the records of interest and Z is the corresponding average of the residuals Z, for each record of interest:

Z, = {E_t - \)R_t - {EPI x R, ) = {E_i - 1 - EPl)R_i The corresponding percentage relative standard error of %EI is then 1 %EI %EIxR ¹ σ_ε (ζ ) is calculated under an assumption that individual data values are uncorrected, and that there is negligible contribution to the standard error from uncertainty in the parameter estimates for the models underlying the estimated enhancement values. So a conservative approach is to double the value of %rse(%ET) when using it to construct confidence intervals for the true value of %EI.

After the step 790, the method 700 concludes (740).

Fig. 10 is a flow diagram illustrating an alternative method 1000 of analysing meteorological observation data obtained from the sites (e.g. 340) within the trial area (e.g. 300) at the specified analysis intervals over the trial period. The method 1000 may be implemented as software resident on the hard disk drive 210 and being controlled in its execution by the processor 205 of the computer system 200.

The method 1000 does not require designation of separate target and control areas within the trial area. This method 1000 is used when a satisfactory control area, free from the influence of the GREA but usefully correlated with the target area, cannot be defined. The method 1000 will be described with reference to a second trial, in South Australia. Some of the main data elements of the second trial are summarised below.

Two GREAs were used in the second trial. The first GREA C2 was located at 35°18' 41.34'S, 138° 31 ' 22.02'E, 44km south-southwest of the Adelaide CBD and approximately seven km from the coast on the Gulf of St Vincent. In addition to this site, a second GREA C3 was located at Tea Tree Gully (34°49' 28.10'S, 138° 44' 48.70'E) around 58km north east of C2. The second GREA C3 was located 18km northeast of the Adelaide CBD and approximately 25km from the coast on the Gulf of St Vincent. The locations of C2 and C3 were sufficiently near each other for there to be a reasonably strong correlation between the natural rainfall in each area of influence. Variations in rainfall were assessed through an analysis of rainfall data from the Bureau of Meteorology rain gauges located within a circle of radius approximately 90 km centred on each GREA site. Although the true extent of any GREA influence is unknown, this trial area was selected on the basis of being likely to capture the bulk of gauges whose measured rainfall, based on previous trials, would be expected to be influenced by the operation of the GREA. There were 282 rain gauges within this trial area. A small number of gauges (20) which did not provide data of sufficient quality were excluded from the trial, leaving 262 gauges that provided rainfall data over the period of the trial. Only 85 of these gauges provided data for every day of the trial.

The GREA sites in the second trial were both in the Mount Lofty Ranges, which are orientated northeast to southwest. Both sites are along the first significant ridgeline closest to the coast, and are exposed to the prevailing weather— typically from the west

South Australia is classified as having a Mediterranean climate and is influenced by offshore trade winds in the summer and on-shore westerlies in the winter. As a consequence, the trial location experiences a dry and warmer period from November to April with prevailing winds from the southeast to east and a moderately wet and colder period from May to October with prevailing winds from the northwest to southwest (BOM, 2008). The climate of the Mount Lofty Ranges is significantly affected by an elevation ranging from 350m to 700m and winds sweeping across the Gulf of St Vincent.

The C2 and C3 sites were located at an elevation of 348m and 373m above sea level respectively, and have significant upslope valleys located to the west and northwest . At C2, the landform elevation rises from the coast travelling from west to east for 4300m at a 1.1 per cent rise (i.e. 1.1m vertical for every 100m horizontal), then continuing east for 2100m the rise increases to 12.3 per cent with the last 200m corresponding to a very steep 21.7 per cent rise. Similarly, C3 has an elevation rise from the coast travelling from west to east for 23,000m at 3.3 per cent followed by a steeper rise of 1 1.2 per cent over the final 2,000m.

Typically, a moist marine onshore airflow from the west rises as the airflow approaches these sites— i.e. there is orographic lifting. The resultant turbulence and vertical movement of air would be expected to result in quick upward dispersal of the ions generated by each GREA.

The second trial ran for 128 days subject to the operating protocol described below, commencing at 9am 1 August 2009 and finishing at 9am 7 December 2009. The analysis interval was 24 hours, from 9 a.m. one day to 9 a.m. the next.

During the second trial, the two GREAs were switched on or off at 9am each day in accordance with a randomised switching regime. This was to coincide with the analysis interval, and to reduce the chance that overlap of rainfall measurements diluted the results. A 30-minute 'temporal buffer' was also added to the switch time, so that the ions from the off-going GREA had time to clear the area before switching on the ongoing GREA. Thus, with a nominal switch time at 9am, the operating GREA was turned off at 8.30am and the ongoing GREA was then turned on at 9am. C2 was operated on a randomised on-off sequence. C3 was operated on a randomised on/on-off/off sequence. To achieve this schedule, consecutive two-day blocks for each site were generated, with 1= on, and 0 = off (see Table 2).

Table 2: Schedule for consecutive two-day blocks.

Auxiliary data for use in the analysis included daily meteorological observations from Adelaide airport and the location and elevation of each rainfall gauge. Observations from Adelaide airport were determined as daily averages and included:

Wind speed (km/h);

Wind direction (degrees from due north, clockwise) with separate readings at 700 hPa, 850 hPa and 925 hPa;

Air temperature;

Dew point temperature; and

Mean sea level pressure.

Steering wind direction and speed for the trial were approximated by a vector average of the 850 and 925 hPa values of wind direction and speed.

Locating an external static control area that matched the trial area, as described above with reference to Fig. 3, was difficult for the second trial. The meteorological and topographic characteristics of neighbouring areas were quite different from those of the trial area. The land area to the north and east of the trial area is relatively flat and dry when compared to the trial area, and the influence of offshore fronts on precipitation in these areas is not nearly as strong.

In the crossover design used for the second trial, the downwind area of C2 acts as a control area for the downwind area of C3 when C2 is off and C3 is on, and vice versa. The "second implementation" described above is used to define downwind areas for each GREA based on steering wind direction in each analysis interval. The effect of "seeding" (GREA operation) can then be assessed using the value of the root double ratio (RDR), which is the geometric mean of the ratios of the area-specific seeded to unseeded rainfall. This can be expressed as:

where

C2_on denotes average rainfall downwind of C2 when C2 is operational and C3 is not;

C2_0ff denotes average rainfall downwind of C2 when C3 is operational and C2 is not;

C3_0n denotes average rainfall downwind of C3 when C3 is operational and C2 is not;

C3_0ff denotes average rainfall downwind of C3 when C2 is operational and C3 is not.

The root double ratio statistic has an expected value of one if seeding has no effect and there is evidence of a positive effect of seeding upon rainfall if its value is significantly greater than one.

A requirement of double-ratio analysis is that the area downwind of C2 does not overlap with the downwind area of C3, and vice versa. However, given the locations of the GREAs and the dynamically defined target and control areas for each GREA as described above, sites may be:

Downwind of C2 and upwind/crosswind of C3; Downwind of C3 and upwind/crosswind of C2; Downwind of C2 and C3; or Upwind/crosswind of C2 and C3. It is therefore desirable in practice to exclude sites downwind of both C2 and C3 in addition to those upwind/crosswind of both C2 and C3. A modified root double ratio - the dynamic double ratio (DDR) - may therefore be determined as:

_{DDR =} D wnwindC2only_C2on/C3off ^ DownwindC3only_C2off/C3o

DownwindC3only_C2on/aoff DownwindC2only_C2off/C3o where DownwindC2onlyc2on/C3off denotes the 'average rainfall' recorded by sites that were downwind of C2 but not of C3 on days when C2 was operational but C3 was not. Similar interpretations hold for the other components of the DDR.

There are a number of ways that the DDR can be calculated, depending on how the concept of 'average rainfall' is defined. Three such definitions are:

The site / analysis interval average of total observed rainfall in all sites in each of the four downwind areas defined in the DDR over the period of the trial;

The simple average of the average rainfall in each downwind area; and The area-weighted average of analysis interval average rainfall in each downwind area, where the area used for each site by analysis interval rainfall reading is the area of the Voronoi polygon centred at the gauge that provided the reading. This is the polygon defined by locations surrounding the gauge that are closer to it than they are to any other gauge in the trial area.

The first definition is, from a statistical perspective, the most efficient but it does give greater weight to days on which there were more gauge-level observations. Further, in determining the accuracy of the estimate, it would be necessary to take into account the spatial correlation between the site observations. The second definition gives equal weight to all analysis intervals, and can be seen as a comparison of estimates of average rainfall in the downwind areas. However, it does not take into account the fact that the spatial distribution of sites is far from uniform. eighting by the area for which a particular gauge is the closest observation, i.e. its Voronoi area, as in the third definition, corrects for the fact that the spatial distribution of sites is not uniform.

In the first step 1010 of the method 1000, the processor 205 uses Restricted Maximum Likelihood (REML) to fit a random effects linear model known as the Observational Model to the rainfall observations in the trial area over the trial period: LogRain_{i t} = a^{1 '} χ. + β^τγ, + λ^Γζ. , + S^rs_{i t} + _Yi + e. , (1) where:

o i indexes sites and t analysis intervals;

o x is a vector of orographic covariates that are specific to site locations;

o y is a vector of meteorological covariates that vary over time;

o z is a vector of covariates related to the operation of the GREA(s);

o s is a vector of dynamically specified site locations;

o γ is a vector of site-specific random effects;

o f is a random error that varies between sites and analysis intervals; and o α, β, λ, and δ are coefficient vectors.

The inclusion of random effects in the Observational Model allows for correlations in the data that are not captured by the model covariates.

The orographic covariates used in the vector x are site elevation and site location (latitude and longitude).

The covariates related to the operation of the GREA(s) making up the vector z include:

• Dummy variables that identify each GREA's operating status during the analysis interval as well as during the previous analysis interval. The dummy variables account for differences in the one day lag structure of the operating schedule for C2 and C3. This is necessary because C2 was operated on a randomly assigned 2-day cycle (on- off), while C3 was operated on a randomly assigned 4-day cycle (on/on-off/off).

• The distance in degrees from a site to each GREA location site relative to the operating status of the GREA during the analysis interval and the previous analysis interval; and • Site angular location (C2<9, C 0) relative to the steering wind direction at each GREA relative to the operating status of the GREA during the analysis interval and the previous analysis interval.

The dynamic specification of site locations in the vector s corresponds to indicators for whether a site is downwind of the GREA(s) during the analysis interval or the previous analysis interval, as well as variables corresponding to the distance(s) of the site from the GREA(s) when the site is downwind, based on the steering wind direction for the analysis interval as determined using the method 800 (see Fig. 8).

The meteorological covariates making up the vector y include, but are not limited to:

• Seasonal variation - represented by a fixed period effect, for example a fixed month effect, a fixed season effect, or a fixed annual effect;.

• Lagged rainfall (expressed as the natural logarithm of observed rainfall in the previous analysis interval);

· Meteorological conditions during the analysis interval - wind direction, wind speed, barometric pressure, and (depending on the length of the trial period) the Southern Oscillation Index (SOI); and

• Meteorological conditions during at least the previous analysis interval, i.e. lagged values of wind direction, wind speed, barometric pressure, air temperature, dew point temperature.

• A fixed indicator variable for Widespread Rainfall Event (WRE) days; and

• A constant.

There were days in the second trial when an exceptional amount of rainfall was recorded throughout the trial area, and not just downwind of one of the GREAs. In particular, on seven of the trial days at least 10mm of rainfall was recorded in at least 250 of the 301 sites that provided data for the trial. These days are defined as WRE days. Such days have a substantial influence on statistical analysis given the squared error loss function of the Observational Model. The inclusion of the WRE indicator removes the average effect of WRE days from the Observational Model. Including a fixed indicator variable for WRE days also improves' the model fit, reducing the variability in the model errors.

Two types of random effects may be included in the random effects vector γ. The first type allocates a random effect to each individual site. This is a spatial random effect that allows for potentially unmeasured orographic variables that vary from one site to the next. The second type of random effect is a spatio-temporal random effect defined by classifying site-level observations during an analysis interval according to their downwind and crosswind orientations relative to the two GREA locations. Inclusion of this spatio- temporal effect reflects the hypothesis that the correlation between site-level rainfall measurements will be stronger in a downwind direction as opposed to a crosswind direction. The six spatio-temporal classifications were defined as described above with reference to Table 1.

The spatio-temporal grouping in accordance with Table 1 leads to 128 (number of analysis intervals of the second trial) x 6 groups = 744 spatio-temporal classifications. However, only (six hundred and sixty four) 664 such groups actually contained site- interval data over the period of the second trial. The large number of spatio-temporal groups and the way in which they were constructed again leads to the possibility that the random effects will attenuate part of the GREA signal. The radial angle and distance to the mid-point of each group may be a proxy for exposure to the GREA ion plume during a given analysis interval. This would lead to an under-attribution if there is in fact a positive GREA contribution to rainfall. In the second trial, only random effects of the second type (spatio-temporal classification) were included. The fit of the Observational Model was 71.1%, indicating that the Observational Model accounts for nearly three-quarters of the variation in LogRain.

Step 1020 follows, in which the processor 205 estimates the enhancement effect of the GREA. The aim of step 1020 is to decompose the observed rainfall for site / on analysis interval t when rainfall is observed at the site as follows:

ObservedRainfall ) , = LatentRainfall _{i t} (l + EnhancementEffect_{i t} ) (2) where LatentRainfall ] , is the natural rainfall that would have been observed at site / in analysis interval t if one or both GREAs had not been operating. Since latent rainfall cannot be observed while the GREAs are operating, the processor 205 determines a prediction of the log scale values of the GREA enhancement effect using the coefficient vectors λ and δ obtained as part of the Observational Model (1), i.e. excluding the variables unrelated to the operation of the GREAs:

LogGREAEffect_{i t} = A^T _{Zi i} + 5^Ts_{i t} (3)

The predicted log scale GREA effects defined by (3) are then mean corrected. This has the effect of moving the expected value of the log scale GREA effects into the corresponding expected value of the log scale latent rainfall.

The processor 205 then determines estimates of the GREA enhancement effect for a particular site-interval when rainfall is observed as

The factor k corrects for the bias that is inherent in using exponentiation to move from log scale rainfall to raw scale rainfall. This bias arises because an effect that changes the mean on the log scale has an asymmetric effect on the variance at the raw scale, understating positive residuals and overstating negative residuals. The method used to calculate the bias correction factor k is now described with reference to Fig. 12.

Fig. 12 is a flow diagram illustrating a method of determining a bias correction factor as used in step 1020 of the method 1000. The method 1200 may be implemented as software resident on the hard disk drive 210 and being controlled in its execution by the processor 205 of the computer system 200. In the first step 1210 of the method 1200, the processor 205 determines a prediction of the logarithm of LatentRainfall_{i t} using the coefficient vectors a, β, and γ obtained as part of the Observational Model (1), by setting

LogLatentRainfall_{i t} = a^Tx_{i t} + β^Τ y_u + At the next step 1220, the processor 205 determines the predicted total rainfall at each site / interval using the predicted values of LogLatentRainfall_{j t} (from step 1210) and

LogGREAEffect_{i t} (from step 1020):

PredictedRainfall_j , =

+ LogGREAEffect_{i t} ) Step 1230 follows, at which the processor 205 determines the variances of the predicted LogLatentRainfall _{i t} and LogGREAEffect_{i t} . Finally, at step 1240, the processor 205 determines the bias correction factor k as

where m is the ratio of the variances determined in step 1230, and r is given by:

ObservedRainfall . ,

r ^■1

n ^ PredictedRainfa.il , , Returning to the method 1000, in step 1030 the processor 205 determines estimates of the contribution of the GREA to rainfall at a site in an analysis interval when rainfall is observed. First, LatentRainfall_{j t} is determined from (2) as ObservedRainfall . ₍

LatentRainfall ^',. , = ^: (5)

EnhancementEffect_i , + 1

The estimated contribution of the GREA to rainfall at a site in an analysis interval when rainfall is observed is then determined as

GREAContribution_{i t} = ObservedRainfall , - LatentRainfall , (6)

In the final step 1040, the processor 205 estimates the precision (standard error) of the total estimated GREA contribution (6) for domains defined by specified site intervals. The calculation of standard errors is based on an assumption of spatial independence of rainfall between sites. Since this assumption is erroneous, the calculated standard errors are inflated by 100 per cent in an attempt to arrive at conservative estimates of the true standard errors. Confidence intervals are then calculated on the basis that errors associated with the estimated GREA contribution are normally distributed. The method 1000 then concludes.

Standard normal theory methods for constructing confidence intervals for complex statistics, as in step 1040, can be sensitive to the assumption of underlying normality. This is especially the case when the complex statistics are based on estimated random effects, as is the case for the method 1000. The distribution of these contribution estimates may be simulated or "bootstrapped" in as non-parametric way as possible in order to get around this problem. A semi-parametric block bootstrap approach may alternatively be used. A block bootstrap simulates the sampling distribution of a complex statistic based on correlated data by resampling blocks of data values rather than individual data values, with the blocks constructed so that, as far as possible, they contain individual values that are correlated with one another within a block, but are uncorrelated across blocks. In the context of the Observational Model, the blocks are the spatio-temporal groups underpinning the random effects, defined in Table 1, and the values associated with each block are its average residual and the within-block deviations from this average that together make up the unconditional residuals for that block under the Observational Model.

In order to construct the semi-parametric block bootstrap distribution for a statistic that depends on the Observational Model, the Observational Model was first used to simulate rain events for all 8661 site-intervals contributing to the results. Working with blocks containing sites designated as recording a rain event following this process, two independent random block samples were selected with replacement, with the first sample contributing the average residuals for these blocks, and the second contributing the within- block deviations for all gauges identified as recording rain in the same blocks. These values were then combined with the estimated fixed effects generated via the Observational Model to produce a set of simulated log rainfall data values for each such block. Finally, these log values were exponentiated in order to _. recover , actual rainfall measurements.

In order to ensure that the rainfall distribution generated in this way is as realistic as possible, three further modifications were made to this distribution - all simulated rainfall values of less than 1 mm were rounded to the closest multiple of 0.2 mm; all simulated daily rainfall values greater than 50 mm were randomly restricted to lie between 25 mm and 50 mm; and the total simulated rainfall for the trial period was randomly restricted to lie between 10,000 mm and 17,000 mm (the actual total for the year 2009 was 13,640 mm). Finally, these simulated rainfall values were used to refit the Observational Model, and all statistics (including contribution estimates) were re-calculated. For the results quoted below, ten thousand (10,000) bootstrap repetitions were used in order to generate bootstrap distributions for the statistics of interest.

The use of a non-parametric block bootstrap, where the blocks correspond to groups in a mixed model, results in bootstrapped estimates of the variance components in this model that are negatively correlated. This leads to substantial undercoverage of then- resulting bootstrap confidence intervals. Accordingly, a post-bootstrap correction was made to the bootstrap distributions obtained using the steps described above. First, the multivariate bootstrap distribution of the variance component estimates for the Observational Model was transformed ('tilted') in order to ensure that these estimates are uncorrected. Next, bootstrap distributions of model parameter estimates were 'tethered' to the original estimate values, using either a mean correction (for estimates, e.g. regression coefficients, defined on the entire real line) or a ratio correction (for estimates, e.g. variance components, that are strictly positive). Finally, bootstrap distributions of complex statistics (including attribution estimates) dependent on these parameters were redetermined. Note that these distributions were also 'tethered' to their original sample values. The contribution estimates defined by equation (6) are unweighted, reflecting only the estimated increase in level of rainfall as measured in the downwind sites that contributed to the modelling process. Since these sites are not distributed uniformly over the trial area, the contribution estimates defined by equation (6) should not be interpreted as estimating the increase in the volume of rain that fell in this area. Accurate conversion of site measurements to volume of rainfall on the ground requires sophisticated spatial modelling and prediction that is beyond the scope of this disclosure. However, an estimate of rainfall volume may be obtained by multiplying each site rainfall observation by the area of the Voronoi polygon surrounding the site and then summing over sites. Note that a Voronoi polygon for a particular site identifies the region containing points that are closer to that site than they are to any other site. Consequently, contributions were calculated on an unweighted and on a Voronoi-area-weighted basis, with the latter providing estimates that are more closely aligned with the volume of rainfall that fell in the trial area. Note, however, that Voronoi weighting tends to give large weights to sites in regions with sparse coverage. As a consequence, the weighted contribution estimates tend to be relatively more variable.

Table 3 summarises the corresponding estimates of latent rainfall based on equation (5) as well as the estimated total GREA contribution based on equation (6) for all site-interval observations to which the Observational Model (1) can be fitted. The overall estimated GREA contribution within the trial area over the trial period is 9.4 per cent.

Table 3: Results summary

The overall standard error of the GREA contribution given in Table 3 was estimated using the semi -parametric block bootstrap method described above was 6.8 per cent.

Table 4 presents bootstrap-based one-sided confidence intervals at different confidence levels for the unweighted and Voronoi weighted contributions (as a percentage of estimated natural rainfall) over the second trial. Table 4 shows that from an unweighted perspective there is ninety (90) per cent confidence that the GREAs made a positive contribution to rainfall over the trial period. This confidence drops to eighty (80) per cent based on Voronoi-weighted or volumetric contribution. This decline is not unexpected given that large weights are given to gauges in sparsely covered areas, and reflects the fact that there is inherently more variability in volumetric rainfall measurements than in site measurements when sites are widely scattered over a large area. As a consequence the precision of any volumetric estimate of rainfall based on Voronoi weighting of site data may not be very high.

Confidence level 99% 95% 90% 80% 70% 60% 50% Unweighted

Atlant attribution (%) -4.8 -1.1 1.0 3.6 5.6 7.4 9.1

Voronoi area weighted

Atlant attribution (%) -8.0 -4.2 -2.3 0.1 2.1 3.7 5.3

Table 4: Lower bounds for parametric bootstrap estimates of one-sided confidence intervals for GREA contribution based on level (unweighted) and volume (V oronoi-area-weighted) Fig. 1 1 is a flow diagram illustrating a method 1 100 of analysing rainfall and other meteorological observations obtained from the sites (e.g. 340) within a trial area, e.g.

300, at specified analysis intervals over a trial period in order to estimate an increase in the probability that a rainfall event occurred during the trial period within the trial area due to the operation of a GREA, e.g. 305. This is the first component of the enhancement effect of the GREA mentioned above. The method 1 100 may be implemented as software resident on the hard disk drive 210 and being controlled in its execution by the processor

205 of the computer system 200.

In the first step 1 1 10 of the method 1100, the processor 205 determines a first model for the probability that a rainfall event did not occur at each site in the control areas for each analysis interval during the trial period using logistic regression. The dependent variable for the first model is a binary variable with a value of one if no rainfall was recorded over the interval, and zero otherwise. The explanatory variables included in the first model are latitude and longitude as well as a fixed day effect. The probability of not observing a rainfall event at each site was adjusted according to the proportion of sites that did not record a rainfall event across the control areas.

In the next step 1120, the processor 205 predicts probabilities of rainfall events for each site in the control and target areas and each analysis interval using the first model. A rainfall event is predicted to not occur at a given site / interval if the first model specifies that the probability of a rainfall event not occurring is greater than 50 per cent.

Step 1130 follows, at which the processor 205 determines a second model for the probability that a rainfall event did not occur at each site in the control and target areas for each analysis interval during the trial period, again using logistic regression. The linear predictor associated with the first model, i.e. the logit of the predicted probability of a rainfall event obtained at step 1 120, is included as an explanatory variable. The second model also includes dummy variables for sites in the north and south control areas, as well as one for sites located in a downwind target area defined by the wind flow sector corresponding to an angle of 30 degrees either side of the downwind direction from the GREA location.

In the final step 1 140, the processor 205 predicts rainfall events for each site in the control and target areas and each analysis interval using the second model. A comparison may be made between the aggregated probability over sites in the target area and that over sites in the control areas to determine whether the operation of the GREA is associated with an increased probability of rainfall events.

Industrial Applicability

The arrangements described are applicable to the agricultural and other weather- sensitive industries.

The foregoing describes only some embodiments of the present invention, and modifications and/or changes can be made thereto without departing from the scope and spirit of the invention, the embodiments being illustrative and not restrictive.

In the context of this specification, the word "comprising" means "including principally but not necessarily solely" or "having" or "including", and not "consisting only of. Variations of the word "comprising", such as "comprise" and "comprises" have correspondingly varied meanings.

Claims

Claims:

1. A method of estimating the effect of a weather modification apparatus over a trial area, the method comprising:

determining an observational model for meteorological observations at sites over the trial area using variables unrelated to the operation of the weather modification apparatus and variables related to the operation of the weather modification apparatus;

determining a first set of meteorological values in the trial area using the observational model excluding the variables unrelated to the operation of the weather modification apparatus; and

determining an estimate of the effect of the weather modification apparatus over the trial area using the first set of meteorological values.

2. A method according to claim 1, wherein the variables unrelated to the operation of the weather modification apparatus include one or more variables selected from the group consisting of:

a fixed period effect;

meteorological observations in the previous analysis interval;

historical meteorological observations;

orographic variables;

meteorological conditions during the analysis interval;

meteorological conditions during the previous analysis interval; and

a random effect.

3. A method according to claim 2 wherein the random effect is a classification of the observation site dependent on the angular location of the observation site relative to the direction downwind of the apparatus during the analysis interval.

4. A method according to claim 1 , wherein the variables related to the operation of the weather modification apparatus include one or more variables selected from the group consisting of:

distance of the observation site from the apparatus;

operation status of the apparatus during the analysis interval and one or more previous analysis intervals;

angular location of the observation site relative to the direction downwind of the apparatus during the analysis interval and one or more previous analysis intervals.

5. A method according to claim 3 or claim 4 wherein the weather modification apparatus emits an ion plume, and the downwind direction is calculated as a speed weighted average of a representative wind profile measured in a layer of the atmosphere considered to contain the majority of winds driving the dispersion of the ion plume.

6. A method according to claim 1, wherein the apparatus is a ground-based rainfall enhancement apparatus and the meteorological observations are rainfall observations.

7. A method according to claim 1, further comprising determining estimates of the contribution of the apparatus to the observed meteorological values over the trial area.

8. A method according to claim 7, further comprising determining the precision of the estimates of the contribution of the apparatus over the trial area.

9. A method according to claim 8, wherein the precision is determined using a semi- parametric block bootstrap method.

10. A computer readable storage medium having a computer program recorded thereon, the program being executable by a computer apparatus to make the computer perform a method of estimating the effect of a weather modification apparatus over a trial area, said program comprising:

code for determining an observational model for meteorological observations at sites over the trial area using variables unrelated to the operation of the weather modification apparatus and variables related to the operation of the weather modification apparatus; code for determining a first set of meteorological values in the trial area using the observational model excluding the variables unrelated to the operation of the weather modification apparatus;

code for determining an estimate of the effect of the weather modification apparatus over the trial area using the first set of meteorological values.

11. A method of estimating the effect of a weather modification apparatus over a trial area comprising at least one control area and a target area, the method comprising:

determining a control model by regressing meteorological observations in the trial area against variables unrelated to the operation of the weather modification apparatus; determining an effects model by regressing the meteorological observations in the trial area against the variables unrelated to the operation of the weather modification apparatus, and variables related to the operation of the weather modification apparatus; determining a first set of meteorological values in the target area using the control model and a second set of meteorological values in the target area using the effects model; and

determining a difference between the first set of meteorological values and the second set of meteorological values, the difference representing an estimate of the effect of the weather modification apparatus over the trial area.

12. A method according to claim 1 1, further comprising:

determining a first regression model of a relationship between meteorological observations in the control areas and locations at which the observations were recorded; and

calculating expected meteorological values over the trial area using the first regression model,

wherein the variables unrelated to the operation of the weather modification apparatus include the expected meteorological values.

13. A method according to claim 1 1, wherein the determining of the first regression model comprises:

determining a trend in the meteorological observations in the control areas;

calculating the deviations of the meteorological observations in the control areas from the estimated trend; and regressing the deviations against the locations in the control areas at which the observations were recorded.

14. A method according to claim 13, wherein the calculating of expected meteorological values comprises adding values obtained from the estimated trend to deviations obtained from the first regression model.

15. A method according to claim 1 1, wherein the variables unrelated to the operation of the weather modification apparatus include one or more variables selected from the group consisting of:

a fixed period effect;

meteorological observations in the previous analysis interval;

historical meteorological observations;

orographic variables;

meteorological conditions during the analysis interval; and

meteorological conditions during at least the previous analysis interval.

16. A method according to claim 1 1, wherein the variables related to the operation of the weather modification apparatus include one or more variables selected from the group consisting of:

activity status of the apparatus;

distance of the observation site from the apparatus;

duration of apparatus operation during the analysis interval and one or more previous analysis intervals; wind flow sector of the observation site, wherein the wind flow sector indicates the angular distance of the observation site from the direction downwind of the apparatus during the analysis interval.

17. A method according to claim 16, wherein the weather modification apparatus emits an ion plume, and the downwind direction is calculated as a speed weighted average of a representative wind profile measured in a layer of the atmosphere considered to contain the majority of winds driving the dispersion of the ion plume.

18. A method according to claim 1 1, wherein the apparatus is a ground-based rainfall enhancement apparatus and the meteorological observations are rainfall observations.

19. A method according to claim 18, further comprising determining an estimate of the volume of rainfall attributable to the apparatus over the trial area.

20. A computer readable storage medium having a computer program recorded thereon, the program being executable by a computer apparatus to make the computer perform a method of estimating the effect of a weather modification apparatus over a trial area comprising at least one control area and a target area, said program comprising:

code for determining a control model by regressing meteorological observations in the trial area against variables unrelated to the operation of the weather modification apparatus;

code for determining an effects model by regressing the meteorological observations in the trial area against the variables unrelated to the operation of the weather modification apparatus, and variables related to the operation of the weather modification apparatus; code for determining a first set of meteorological values in the target area using the control model and a second set of meteorological values in the target area using the effects model; and

code for determining a difference between the first set of meteorological values and the second set of meteorological values, the difference representing an estimate of the effect of the weather modification apparatus over the trial area.