Executive Summary

This report details the process of calibrating and validating a water balance model for the Melbourne catchment, to obtain daily soil moisture values and runoff. The model uses a split sample, whereby BoM precipitation data values from 2005-2009 were used to calibrate the model and values from 2011-2016 were used to validate the model. The simple water balance considers the daily change in storage that is a function of precipitation, evapotranspiration, surface runoff and baseflow. The model was constructed with the following simplifying assumptions; soil depth is 785 mm, wilting point is 0.12m³m^-3, field capacity is 0.33 m³m^-3, and saturation is 0.45 m³m^-3. Actual evapotranspiration was considered in the model, with daily areal potential evapotranspiration values obtained from BoM.

The model was calibrated to find a suitable value for runoff coefficient (RC) and k (baseflow factor), by linearly spacing 100 values between 0-1 for RC and 0-1 for k, running the model and testing sensitivity through the Nash-Sutcliffe Efficiency (NSE) Coefficient. From the calibration, a value of 0.0182 for k and 0.0707 for RC were determined. The water balance model was then updated to reflect the new constant values for k and RC and implemented for data from 2010-2016. Results of the simulation were compared against known runoff and soil moisture data from BoM for 2011-2016 to validate the accuracy of the model with the calibrated parameters.

After comparing the water balance model outputs with the known values from BoM, it was determined that the calibrated parameters were valid when comparing with the runoff from 2011-2016, as the model roughly followed similar patterns to the actual runoff measured. When comparing the soil moisture values, however, the observed and simulated values differed significantly, which is likely attributed to the simplifying assumptions used to initially create the model, such as a constant soil depth, wilting point and maximum saturation. Therefore, there can be a higher level of confidence in using the model to simulate runoff in the Melbourne catchment than exists for storage and soil moisture values.

Executive Summary i

Introduction 2

1. Water Balance Model 2

1.1. Finding ETa 2

1.2. Finding Total Runoff 3

2. Water Balance Model Calibration 3

2.1. Calibration Process 3

2.2. Calibration Result 4

3. Water Balance Model Validation 8

3.1. Validation Results 8

Conclusion 10

References 10

Appendix A – Calibration & Validation Code 10

Introduction

Water management is one of the important parts of the water cycle. This means to plan, develop, distribute and optimise the necessary usage of the total water resource especially for modern society with a rising population around the world. Poor water management can lead to a huge waste of water, and serious greenhouse effects that can have adverse effects on the environment, leading to climate change. In the future, predicting the water availability will be a key step goal for maintaining Biospheric life, including livestock, humans and the entire ecosystem consumption.

As a semi-arid country, the Australian government set up a system called Bureau of Meteorology (BoM) which is responsible for providing weather forecast and observation services for Australia and surrounding places. Runoff prediction is also one of the functions of the BoM. In this project, our team develops a water balance model for Melbourne with the given information of daily precipitation (P) and evapotranspiration from the BoM website to estimate the runoff for 11 years during 2005-2016 periods. The calibration of k was required using a simple calibration approach as the baseflow was calculated as the product of k and the storage of water. This task however required the calculation of daily evapotranspiration values by the water balance model achieved by finding ET_a which was done by utilising the relationship between ET_p and S. Various assumptions were made throughout the task including a reasonable value for the storage of water in soil based on location and starting time of the simulation. The soil depth was assumed to be 785m, wilting point at 0.12, saturation 0.4 and field capacity at 0.33 (all values in units of m^3 m^-3).

Regarding the wilting point, it was believed that once the wilting point was achieved, water cannot be further depleted from the water balance model and once the storage was full, no further water can enter the saturated soil. Parameters in the water balance models such as field capacity, total size and wilting point were tested by individually analysing their response to the model. Moreover, the task regarding calibration and validation was influential to the validity of this task and the standard split-sample test was used which allowed for the utilisation of 0.33-0.5 of the data for calibration and the others for validation. The Monte-Carlo simulation approach was used where the model was run at least 100 times with each combination of parameters that were used because of this. The Nash-Sutcliffe coefficient was then stored in a 100×100 matrix each time thus calculating the objective function which was then run with the best value amongst the whole set of runs. After receiving the Bureau of Meteorology’s soil moisture data for the same period, the found runoff and storage terms were compared for the two methods. With this information, the model was thoroughly analysed for validity and accuracy whilst creating different runs.

Water Balance Model

The water balance model in this report is based on Equation 1, where P is precipitation (mm/day), ETa is evapotranspiration (mm/day), Rt is the total runoff (mm/day).

Precipitation values are known from rainfall observation data, however Eta and Rt are unknown.

Finding ETa

To find the contributing daily value of actual evapotranspiration (ETa), potential evapotranspiration (ETp) data was sourced from the Australian Bureau of Meteorology. ETp can be defined as the potential amount of evapotranspiration that may occur, assuming that conditions are neither energy nor water limited. In practice, it is unlikely that both conditions will be reached, thus the ETa must be determined, taking into consideration the limiting factors for each day. ETa can be found from Equation 1.1 .1.

Equation 1.1.1

Thus, to find ETa, it is required to first find the value of SI, by using Equation 1.1 .2.

Equation 1.1.2

The field capacity is known to be 0.33m³m^-3 and the wilting point is known to be 0.12 m³m^-3. Soil moisture can be found using Equation 1.1 .3.

Equation 1.1.3

Soil depth is known to be 0.785m, however the daily storage values remain unknown and will be found through running the model.

Finding Total Runoff

The total daily runoff in the catchment consists of both surface runoff and baseflow contributions. The total runoff can be defined by Equation 1.2 .4, where Rs is surface runoff and Rb is baseflow.

Equation 1.2.4

Surface runoff (Rs) can be found through Equation 1.2 .5, where RC is runoff coefficient and P is precipitation.

Equation 1.2.5

The runoff coefficient, RC, must be calibrated using the water balance model.

Baseflow (Rb) can be found through Equation 1.2 .6.

Equation 1.2.6

Daily storage data, S, can be found through running the water balance model, however the model must be calibrated to find k.

The total runoff in the catchment increases beyond these values when the soil is at saturation and there is further excess runoff. Excess runoff has been calculated by Equation 1.2.4, when storage is higher than the maximum storage.

Equation 1.2.7

Water Balance Model Calibration

Calibration Process

The water balance model must be calibrated such that a constant value for both RC and k are determined. Once calibrated, the model can be used to produce daily change in storage data for the Melbourne catchment, including the daily runoff and actual evapotranspiration. To calibrate the model, it must first be initialised with an arbitrary storage value. The initial storage value chosen was the midpoint of the wilting point storage and maximum saturation storage, as this represents the midrange value for storage in the model.

The model has been limited based on maximum and minimum storage. When storage is higher than the maximum storage of 0.45m3/m3, no more water can enter the soil matrix. Thus, excess water is converted to runoff and added to the daily runoff value. The minimum storage of the model is set to 0, thus if storage of a particular day is negative, it will be set to 0. For evapotranspiration to occur, the storage must be above the wilting point. Thus, the daily evapotranspiration is set to 0 if the storage is less than or equal to the wilting point.

The Monte-Carlo simulation approach was employed in the model, whereby 100 evenly spaced values of k from 0 to 1 and 100 evenly spaced values of RC from 0.09-0.71 were trialled against the objective function, with each combination of both considered. The RC range of 0.09-0.71 was chosen as Melbourne has a temperate climate with wet winters, classified as Cs, which has a 10^th percentile runoff coefficient of 0.09 and a 90^th percentile coefficient of 0.71. However, upon trialling the simulation, it was found that the RC given from the approach was 0.09, the absolute lower bound of the range. Thus, the range was expanded from 0-1 to ensure that the full range of possible RC values were tested.

The model then utilised the Nash-Sutcliffe Efficiency (NSE) Coefficient as the objective function to determine the sensitivity of the model for different values of k and RC. The NSE compares given runoff values from BoM and the runoff values obtained from the model, for each combination of k and RC, using Equation 2.1 .8.

Equation 2.1.8

Calibration Result

The results of the NSE were plotted against each k and RC value to produce a surface, as seen in Figure 2.2 .1, and contour plot, as seen in Figure 2.2 .2, to determine an appropriate value of k and RC for the model.

Figure 2.2.1: NSE Value Surface Plot for Full Range of k and RC

Figure 2.2.2: Contour Graph of NSE Values for Full Range of k and RC

Using Figure 2.2 .2 as a reference, the ranges of k and RC were adjusted to the bounds of the maximum yellow contour, such that k was assigned a range of 0-0.1 and RC assigned 0-0.2. By reducing the ranges of the parameters, a higher resolution can be achieved for a more accurate k and RC value. The resultant Surface plot and contour plot can be seen in Figure 2.2 .3 and Figure 2.2 .4, respectively.

Figure 2.2.3: NSE Surface Plot for Refined k and RC Range

Figure 2.2.4: NSE Contour Plot for Refined k and RC Range

Table 2.2 .1 details the interpretations of different NSE values from 0-1. If the NSE value were smaller than 0, this would indicate that the observed value was a better predictor than the simulated value thus unacceptable performance. The maximum NSE value obtained from the surface in Figure 2.2 .3: NSE Surface Plot for Refined k and RC Range was found to be 0.4234, which lies between the range of 0.36 and 0.75, thus can be interpreted as a “qualified” prediction. While this does not fall within the desired range of 0.75-1, the calibrated parameter combination can be proceeded with.

Table 2.2.1: NSE Value Interpretation [3]

NSE Value	Interpretation
0.75	Good
0.36	Qualified
NSE	Not Qualified

The optimum runoff coefficient and k value that correspond to the maximum NSE value are 0.0707 and 0.0182, respectively. While the calibration process assumes that performing the model with these parameters for the 2005-2009 data will provide a similar result for runoff and soil moisture to the observed values from BoM data, these data sets were plotted in Figure 2.2 .5 and Figure 2.2 .6 to identify whether there are any issues in the model set-up or the simplifying assumptions that may affect the result.

Figure 2.2.5: Calculated Runoff v. Observed Runoff from BoM

Figure 2.2.6: Calculated Soil Moisture v. Observed Soil Moisture from BoM

From these results, it can be determined that the runoff calculation has provided a successful result that roughly follows the observed runoff, which is to be expected as NSE was calculated based on the comparison of calculated and observed runoff. It should be noted that runoff has been overestimated by the model for some time periods, which may suggest that excess runoff is overestimated due to an underestimation in maximum storage. The soil moisture has yielded calculated values that while following the pattern of observed values, are approximately 0.6-0.7 lower, hence providing a large margin of error. This difference may be attributed to incorrect simplifying assumptions on the catchments soil depth, wilting point, field capacity and maximum storage. Thus, further research should be performed to obtain values that provide a better representation of the Melbourne catchment.

Water Balance Model Validation
1. Validation Results

The calibrated parameters from the above section were tested against an independent data set to determine whether the parameters found were accurate in producing a water balance model for the Melbourne catchment beyond the initial calibration range of 2005-2009. The water balance model was applied to precipitation data from the Melbourne catchment for 2010-2016 and then compared with data provided from the Bureau of Meteorology for runoff and soil moisture. Figure 1.5.1 shows a plot of Observed Runoff from BoM and Calculated Runoff from the simulation for 2011-2016. The graph shows that the simulation follows the same peaks and troughs as the observed runoff through the period. The calculated runoff does not experience the same extreme highs as the observed runoff. This may suggest that the way in which excess runoff has been calculated in the simulation may be underestimating the actual runoff occurring. This may have occurred if the maximum saturation of the catchment was overestimated, thus the excess runoff would be underestimated.

Figure 3.1.7: Calculated Runoff v. Observed Runoff from BoM for 2010-2016

Figure 3.1.8: Calculated Soil Moisture v. Observed Soil Moisture from BoM for 2010-2016

Discussion of results

FOR BOM ETP DATA

Conclusion

From this report, it can be concluded that the calibrated parameters of k and RC for the Melbourne catchment provide an invalid representation of the runoff and soil moisture when applied to 2011-2016 precipitation and potential evapotranspiration data. Thus, it has been determined that the presented method of calibration in this report must be further refined to provide parameter estimations that better model the hydrological climate of the Melbourne catchment.

References

[1] Bureau of Metereology. 2021. Areal potential evapotranspiration, 2005. [online] Available at: [Accessed 31 May 2021].

[2] Bureau of Meterology. 2021. Root zone soil moisture, 1 January 2005. [online] Available at: [Accessed 31 May 2021].

[3] Research Gate. 2021. Criteria of Nash-Sutcliffe Efficiency (NSE) Value. [online] Available at: [Accessed 31 May 2021].

Appendix A – Calibration & Validation Code

clear all; close all; clc;

load(‘qtot_full_A2_MEL.mat’);

load(‘rain_A2_MEL.mat’);

GivenPrec = BoM_other(1:1826);

ValPrec=BoM_other(1827:4383);

Givenrunoff = BoM_runoff(1:1826);

Valrunoff=BoM_runoff(1827:4383);

load(‘sm_full_A2_MEL.mat’);

GivenSm = BoM_other(1:1826);

GivenSMvalid = BoM_other(1827:4383);

% calibration ETp data

data05 = readtable(‘ma_wet_-37.8_145_2005_Actual.csv’);

data06 = readtable(‘ma_wet_-37.8_145_2006_Actual.csv’);

data07 = readtable(‘ma_wet_-37.8_145_2007_Actual.csv’);

data08 = readtable(‘ma_wet_-37.8_145_2008_Actual.csv’);

data09 = readtable(‘ma_wet_-37.8_145_2009_Actual.csv’);

% validation ETp data

data010 = readtable(‘ma_wet_-37.8_145_2010_Actual.csv’);

data011 = readtable(‘ma_wet_-37.8_145_2011_Actual.csv’);

data012 = readtable(‘ma_wet_-37.8_145_2012_Actual.csv’);

data013 = readtable(‘ma_wet_-37.8_145_2013_Actual.csv’);

data014 = readtable(‘ma_wet_-37.8_145_2014_Actual.csv’);

data015 = readtable(‘ma_wet_-37.8_145_2015_Actual.csv’);

data016 = readtable(‘ma_wet_-37.8_145_2016_Actual.csv’);

ObsDates = dates(1:1826); % calibration dates

fpath = ‘C:UserslaureOneDriveDocumentsCIV3285Assignment 2Figures’;

% calibration ETp

ETpdaily = table2array([data05(:,2);data06(:,2);data07(:,2);data08(:,2);data09(:,2)]);

%Given parameters m^3 / m^3

wp = 0.12;

fc = 0.33;

depth = 785;

sat = 0.45;

Swp= wp*depth;

Satdepth = sat*depth;

Depthdiff = Satdepth – Swp;

%Initializing Matrix

S = zeros(100,100,size(Givenrunoff,1)); %Storage

SI = zeros(size(Givenrunoff,1),1);

S(:,:,1) = (Satdepth+Swp)/2; %mm

Runoff = zeros(100,100,size(Givenrunoff,1));

Excess = zeros(100,100,size(Givenrunoff,1));

ETa = zeros(100,100,size(Givenrunoff,1));

dS = zeros(100,100,size(Givenrunoff,1));

NSE = zeros(100,100);

Runoff(:,:,1) = 0;

Excess(:,:,1) = 0;

rangeofrc = linspace(0,0.21,100); %Runoff coefficient – initially run for linspace(0,1,100)

rangeofk = linspace(0,0.17,100); % k value – initially run for linspace(0,1,100)

for day=1:size(Givenrunoff,1)

for k=1:100

for rc=1:100

Runoff(k,rc,day)=rangeofk(rc)*S(k,rc,day)+rangeofrc(k)*GivenPrec(day)+Excess(k,rc,day);

SI(k,rc,day)=S(k,rc,day)/(fc-wp);

ETa(k,rc,day) = min(ETpdaily(day),SI(k,rc,day)*ETpdaily(day));

dS(k,rc,day) = GivenPrec(day) – ETa(k,rc,day) – Runoff(k,rc,day);

S(k,rc,day+1)=S(k,rc,day)+dS(k,rc,day);

if S(k,rc,day+1)

S(k,rc,day+1)= 0;

elseif S(k,rc,day+1)> Satdepth % Making sure the Storage is not higher than the maximum

Excess(k,rc,day+1)= S(k,rc,day+1) – (Satdepth);

S(k,rc,day+1)= Satdepth;

end

% For Objective Function using NSE

for k=1:100

for rc=1:100

Requiredrunoff = squeeze(Runoff(k,rc,1:size(Givenrunoff,1)));

Averagerunoff = mean(Givenrunoff(1:size(Givenrunoff,1)));

ReqGivenrunoff = Givenrunoff(1:size(Givenrunoff,1));

NSE(k,rc)=1-sum((Requiredrunoff-ReqGivenrunoff).^2) / sum((ReqGivenrunoff-Averagerunoff).^2);

end

%Finding the best NSE with the corresponding K and rc values

[rcindex,kindex] = find(NSE==max(max(NSE)));

Koptimum = rangeofk(kindex);

RCoptimum = rangeofrc(rcindex);

S_validation(1)=S(rcindex,kindex,size(Givenrunoff,1));

S_test(1)=S(rcindex,kindex,size(Givenrunoff,1));

% 3D Surface Plot

surf(rangeofk,rangeofrc,NSE)

xlabel(‘k’)

ylabel(‘rc’)

zlabel(‘NSE values’)

title(‘NSE values based on varying K and runoff coefficients’)

print(‘-dpng’,’-r200′,fullfile(fpath,’NSE Restricted’));

% 2D Contour Graph

figure

contour(rangeofk,rangeofrc,NSE)

xlabel(‘k’)

ylabel(‘rc’)

title(‘Contours of varying K and runoff coefficients to NSE – Full Range k and RC’)

print(‘-dpng’,’-r200′,fullfile(fpath,’Contours Restricted’));

%% Ensure calibration works with own data

Runoff_test = zeros(size(Givenrunoff,1),1);

Excess_test = zeros(size(Givenrunoff,1),1);

SMcalc_test = zeros(size(Givenrunoff,1),1);

ETa_test = zeros(size(Givenrunoff,1));

dS_test = zeros(size(Givenrunoff,1));

for day = 1:size(Givenrunoff,1)

Runoff_test(day)=Koptimum*S_test(day)+RCoptimum*GivenPrec(day)+Excess_test(day);

SI(day)=S_test(day)/(fc-wp);

ETa_test(day) = min(ETpdaily(day),SI(day)*ETpdaily(day));

dS_test(day) = GivenPrec(day) – ETa_test(day) – Runoff_test(day);

S_test(day+1)=S_test(day)+dS_test(day);

if S_test(day+1)

S_test(day+1)= 0;

elseif S_test(day+1)> Satdepth %Checking if the Storage …

Excess_test(day+1)= S_test(day+1) – (Satdepth);

S_test(day+1)= Satdepth;

end

SMcalc_test(day) = S_test(day)/depth;

end

figure

plot(ObsDates,Givenrunoff)

hold on

plot(ObsDates,Runoff_test)

xlabel(‘Dates’)

ylabel(‘Runoff (mm)’)

legend(‘Observed Runoff’,’Calculated Runoff’)

title(‘Calibration Runoff 2005-2009’)

print(‘-dpng’,’-r200′,fullfile(fpath,’Runoff Calibration’));

%Storage Calibration

figure

plot(ObsDates,SMcalc_test)

hold on

plot(ObsDates,GivenSm)

xlabel(‘Dates’)

ylabel(‘Soil Moisture (mm/mm)’)

title(‘Storage Calibration 2005-2009’)

legend(‘Calculated Storage’,’Observed Storage’)

print(‘-dpng’,’-r200′,fullfile(fpath,’Calibration Storage’))

%% Validation

%Initializing Matrixes for the Validation Part

Runoff_valid = zeros(size(Valrunoff,1),1);

Excess_valid = zeros(size(Valrunoff,1),1);

SMcalc_valid = zeros(size(Valrunoff,1),1);

ETa_validation = zeros(size(Valrunoff,1));

dS_validation = zeros(size(Valrunoff,1));

ValidationDates = dates(1827:4383);

ETpVal = table2array([data010(:,2);data011(:,2);data012(:,2);data013(:,2);data014(:,2);data015(:,2);data016(:,2)]);

% Validation

for day = 1:size(Valrunoff,1)

Runoff_valid(day)=Koptimum*S_validation(day)+RCoptimum*ValPrec(day)+Excess_valid(day);

SI(day)=S_validation(day)/(fc-wp);

ETa_validation(day) = min(ETpVal(day),SI(day)*ETpVal(day));

dS_validation(day) = ValPrec(day) – ETa_validation(day) – Runoff_valid(day);

S_validation(day+1)=S_validation(day)+dS_validation(day);

if S_validation(day+1)

S_validation(day+1)= 0;

elseif S_validation(day+1)> Satdepth %Checking if the Storage …

Excess_validation(day+1)= S_validation(day+1) – (Satdepth);

S_validation(day+1)= Satdepth;

end

SMcalc_valid(day) = S_validation(day)/depth;

end

%Runoff validation

figure

plot(ValidationDates,Valrunoff)

hold on

plot(ValidationDates,Runoff_valid)

xlabel(‘Dates’)

ylabel(‘Runoff (mm)’)

legend(‘Observed Runoff’,’Calculated Runoff’)

title(‘Validation Runoff 2010-2016’)

print(‘-dpng’,’-r200′,fullfile(fpath,’Runoff’));

%Storage Validation

figure

plot(ValidationDates,SMcalc_valid)

hold on

plot(ValidationDates,GivenSMvalid)

xlabel(‘Dates’)

ylabel(‘Soil Moisture (mm/mm)’)

title(‘Validation Soil Moisture 2010-2016’)

legend(‘Calculated Storage’,’Observed Storage’)

print(‘-dpng’,’-r200′,fullfile(fpath,’Storage Validation’));

Executive Summary

Contents

Introduction

Conclusion

References

Appendix A – Calibration & Validation Code