Merge pull request #14 from TOSSHtoolbox/ry-patch

Streamlined RecessionAnalysis signature output and added more description

Author: RY4GIT

TOSSH_code/calculation_functions/calc_All.m

@@ -106,7 +106,9 @@
 Q_var_error_str = strings(size(Q_mat,1),1);
 QP_elasticity = NaN(size(Q_mat,1),1);
 QP_elasticity_error_str = strings(size(Q_mat,1),1);
-RecessionParameters = NaN(size(Q_mat,1),2);
+RecessionParameters_a = NaN(size(Q_mat,1),1);
+RecessionParameters_b = NaN(size(Q_mat,1),1);
+RecessionParameters_T0 = NaN(size(Q_mat,1),1);
 RecessionParameters_error_str = strings(size(Q_mat,1),1);
 RecessionK_early = NaN(size(Q_mat,1),1);
 RecessionK_early_error_str = strings(size(Q_mat,1),1);
@@ -196,8 +198,17 @@
     [Q_skew(i),~,Q_skew_error_str(i)] = sig_Q_skew(Q_mat{i},t_mat{i});
     [Q_var(i),~,Q_var_error_str(i)] = sig_Q_var(Q_mat{i},t_mat{i});
     [QP_elasticity(i),~,QP_elasticity_error_str(i)] = sig_QP_elasticity(Q_mat{i},t_mat{i},P_mat{i});
-    [RecessionParameters(i,:),~,~,RecessionParameters_error_str(i)] = ...
+    [RecessionParametersTemp,~,~,RecessionParameters_error_str(i)] = ...
         sig_RecessionAnalysis(Q_mat{i},t_mat{i},'fit_individual',false);
+
+    % Post-processing of RecessionParameters (see sig_RecessionAnalysis
+    % function description for the details)
+    RecessionParameters_a(i) = median((RecessionParametersTemp(:,1)),'omitnan');
+    RecessionParameters_b(i) = median(RecessionParametersTemp(:,2),'omitnan');
+    RecessionParametersT0Temp = 1./(RecessionParametersTemp(:,1).*median(Q_mat{i}(Q_mat{i}>0),'omitnan').^(RecessionParametersTemp(:,2)-1));
+    ReasonableT0 = and(RecessionParametersTemp(:,2)>0.5,RecessionParametersTemp(:,2)<5);
+    RecessionParameters_T0(i) = median(RecessionParametersT0Temp(ReasonableT0),'omitnan');
+
     [RecessionK_early(i),~,RecessionK_early_error_str(i)] = sig_RecessionParts(Q_mat{i},t_mat{i},'early');
     [Spearmans_rho(i),~,Spearmans_rho_error_str(i)] = sig_RecessionUniqueness(Q_mat{i},t_mat{i});
     [ResponseTime(i),~,ResponseTime_error_str(i)] = sig_ResponseTime(Q_mat{i},t_mat{i},P_mat{i});
@@ -278,7 +289,9 @@
 results.Q_var_error_str = Q_var_error_str;
 results.QP_elasticity = QP_elasticity;
 results.QP_elasticity_error_str = QP_elasticity_error_str;
-results.RecessionParameters = RecessionParameters;
+results.RecessionParameters_a = RecessionParameters_a;
+results.RecessionParameters_b = RecessionParameters_b;
+results.RecessionParameters_T0 = RecessionParameters_T0;
 results.RecessionParameters_error_str = RecessionParameters_error_str;
 results.RecessionK_early = RecessionK_early;
 results.RecessionK_early_error_str = RecessionK_early_error_str;

TOSSH_code/calculation_functions/calc_McMillan_Groundwater.m

@@ -1,9 +1,9 @@
 function [results] = calc_McMillan_Groundwater(Q_mat, t_mat, P_mat, PET_mat, varargin)
 %calc_McMillan_Groundwater calculates various groundwater signatures.
-%   Calculates 15 signatures from McMillan (2020), related to groundwater 
-%   storage, groundwater dynamics and baseflow. These signatures come from 
-%   previous experimental studies that link catchment or hillslope 
-%   processes to streamflow response dynamics. Some signatures are 
+%   Calculates 15 signatures from McMillan (2020), related to groundwater
+%   storage, groundwater dynamics and baseflow. These signatures come from
+%   previous experimental studies that link catchment or hillslope
+%   processes to streamflow response dynamics. Some signatures are
 %   implemented direct from the original papers, others are interpreted
 %   from a qualitative description in the paper.
 %
@@ -13,7 +13,7 @@
 %   P_mat: precipitation [mm/timestep] matrix (cell array)
 %   PET_mat: pot. evapotranspiration [mm/timestep] matrix (cell array)
 %   OPTIONAL
-%   start_water_year: first month of water year, default = 10 (October)  
+%   start_water_year: first month of water year, default = 10 (October)
 %   plot_results: whether to plot results, default = false
 %   recession_length: min. length of recessions [days], default = 15
 %   n_start: days to be removed after start of recession
@@ -75,7 +75,7 @@
 
 % initialise arrays
 
-% Section: Groundwater 
+% Section: Groundwater
 % Signatures referring to double peaks in streamflow are not coded due to
 % subjectivity in identification of the peaks.
 % Total runoff ratio describes overall water loss to deep groundwater.
@@ -98,73 +98,79 @@
 % Seasonal variations in recession parameters (Shaw and Riha, 2012).
 Recession_a_Seasonality = NaN(size(Q_mat,1),1);
 Recession_a_Seasonality_error_str = strings(size(Q_mat,1),1);
-% Average storage from average baseflow and storage-discharge relationship 
+% Average storage from average baseflow and storage-discharge relationship
 % (McNamara et al., 2011).
 AverageStorage = NaN(size(Q_mat,1),1);
 AverageStorage_error_str = strings(size(Q_mat,1),1);
 % Recession analysis parameters approximate storage-discharge relationship.
 % This fits a single relationship to the point cloud.
-RecessionParameters = NaN(size(Q_mat,1),3);
+RecessionParameters_a = NaN(size(Q_mat,1),1);
+RecessionParameters_b = NaN(size(Q_mat,1),1);
+RecessionParameters_T0 = NaN(size(Q_mat,1),1);
 RecessionParameters_error_str = strings(size(Q_mat,1),1);
-% Changes of slope in a master recession curve (MRC) or recession analysis 
+% Changes of slope in a master recession curve (MRC) or recession analysis
 % plot indicate multiple linear reservoirs.
 MRC_num_segments = NaN(size(Q_mat,1),1);
 MRC_num_segments_error_str = strings(size(Q_mat,1),1);
 % First: steep section of the master recession curve = storage that is
 % quickly depleted (Estrany et al., 2010).
 First_Recession_Slope = NaN(size(Q_mat,1),1);
-% Second: mid section of the master recession curve = water retention  
+% Second: mid section of the master recession curve = water retention
 % capacity of the catchment (Estrany et al., 2010).
 Mid_Recession_Slope = NaN(size(Q_mat,1),1);
 % Non-uniqueness in the storage-discharge relationship (McMillan et al.,
-% 2011; Harman et al., 2009). Tested using Spearman rank correlation on Q 
+% 2011; Harman et al., 2009). Tested using Spearman rank correlation on Q
 % vs. dQ/dt.
 Spearmans_rho = NaN(size(Q_mat,1),1);
 Spearmans_rho_error_str = strings(size(Q_mat,1),1);
-% Ratio between event and total runoff coefficient: low = large storage 
+% Ratio between event and total runoff coefficient: low = large storage
 % capacity (Blume et al., 2008).
 EventRR_TotalRR_ratio = NaN(size(Q_mat,1),1);
-% Variability index of flow. Low variability index shows higher water 
+% Variability index of flow. Low variability index shows higher water
 % storage (Estrany et al., 2010).
 VariabilityIndex = NaN(size(Q_mat,1),1);
 VariabilityIndex_error_str = strings(size(Q_mat,1),1);
 
-% Section: Baseflow 
+% Section: Baseflow
 % Visual inspection of hydrograph for stable baseflow: not coded.
-% Baseflow index indicates baseflow proportion and baseflow residence time 
+% Baseflow index indicates baseflow proportion and baseflow residence time
 % (Yilmaz et al., 2008; Bulygina et al., 2009; Hrachowitz et al., 2014).
 BFI = NaN(size(Q_mat,1),1);
 BFI_error_str = strings(size(Q_mat,1),1);
 % Baseflow recession constant K (assuming exponential recession behaviour),
-% slower recessions show greater groundwater influence and longer 
+% slower recessions show greater groundwater influence and longer
 % subsurface flow paths (Safeeq et al., 2013).
 BaseflowRecessionK = NaN(size(Q_mat,1),1);
 BaseflowRecessionK_error_str = strings(size(Q_mat,1),1);
 
 % loop over all catchments
 for i = 1:size(Q_mat,1)
-    
+
     % Section: Groundwater
     [TotalRR(i),~,TotalRR_error_str(i)] = sig_TotalRR(Q_mat{i},t_mat{i},P_mat{i});
     [RR_Seasonality(i),~,RR_Seasonality_error_str(i)] = sig_RR_Seasonality(Q_mat{i}, t_mat{i}, P_mat{i}, ...
         'summer_start', start_water_year-6);
     [EventRR(i),~,EventRR_error_str(i)] = sig_EventRR(Q_mat{i},t_mat{i},P_mat{i});
     [StorageFraction(i,1),StorageFraction(i,2),StorageFraction(i,3),~,StorageFraction_error_str(i)] = ...
         sig_StorageFraction(Q_mat{i},t_mat{i},P_mat{i},PET_mat{i},'plot_results',plot_results);
-    
+
     % Section: Storage (especially groundwater)
     [Recession_a_Seasonality(i),~,Recession_a_Seasonality_error_str(i)] = ...
         sig_SeasonalVarRecessions(Q_mat{i},t_mat{i},'eps',eps,'recession_length',recession_length,'plot_results',plot_results,'n_start',n_start);
     [AverageStorage(i),~,AverageStorage_error_str(i)] = ...
         sig_StorageFromBaseflow(Q_mat{i},t_mat{i},P_mat{i},PET_mat{i},'start_water_year',start_water_year,'plot_results',plot_results,'recession_length',recession_length,'n_start',n_start,'eps',eps);
     [RecessionParametersTemp,~,~,RecessionParameters_error_str_temp] = ...
         sig_RecessionAnalysis(Q_mat{i},t_mat{i},'fit_individual',true,'plot_results',plot_results,'recession_length',recession_length,'n_start',n_start,'eps',eps);
-    RecessionParameters(i,1) = median((RecessionParametersTemp(:,1)),'omitnan');
-    RecessionParameters(i,2) = median(RecessionParametersTemp(:,2),'omitnan');  
+
+    % Post-processing of RecessionParameters (see sig_RecessionAnalysis
+    % function description for the details)
+    RecessionParameters_a(i) = median((RecessionParametersTemp(:,1)),'omitnan');
+    RecessionParameters_b(i) = median(RecessionParametersTemp(:,2),'omitnan');
     RecessionParametersT0Temp = 1./(RecessionParametersTemp(:,1).*median(Q_mat{i}(Q_mat{i}>0),'omitnan').^(RecessionParametersTemp(:,2)-1));
     ReasonableT0 = and(RecessionParametersTemp(:,2)>0.5,RecessionParametersTemp(:,2)<5);
-    RecessionParameters(i,3) = median(RecessionParametersT0Temp(ReasonableT0),'omitnan');
+    RecessionParameters_T0(i) = median(RecessionParametersT0Temp(ReasonableT0),'omitnan');
     RecessionParameters_error_str(i) = RecessionParameters_error_str_temp;
+
     [MRC_num_segments(i),Segment_slopes,~,MRC_num_segments_error_str(i)] = ...
         sig_MRC_SlopeChanges(Q_mat{i},t_mat{i},'plot_results',plot_results,'eps',eps,'recession_length',recession_length,'n_start',n_start);
     First_Recession_Slope(i) = Segment_slopes(1);
@@ -174,12 +180,12 @@
     [Spearmans_rho(i),~,Spearmans_rho_error_str(i)] = sig_RecessionUniqueness(Q_mat{i},t_mat{i},'eps',eps,'recession_length',recession_length,'n_start',n_start);
     EventRR_TotalRR_ratio(i) = EventRR(i)/TotalRR(i);
     [VariabilityIndex(i),~,VariabilityIndex_error_str(i)] = sig_VariabilityIndex(Q_mat{i},t_mat{i});
-    
+
     % Section: Baseflow
     [BFI(i),~,BFI_error_str(i)] = sig_BFI(Q_mat{i},t_mat{i},'method','UKIH');
     [BaseflowRecessionK(i),~,BaseflowRecessionK_error_str(i)] = ...
-        sig_BaseflowRecessionK(Q_mat{i},t_mat{i},'eps',eps,'recession_length',recession_length,'n_start',n_start); 
-    
+        sig_BaseflowRecessionK(Q_mat{i},t_mat{i},'eps',eps,'recession_length',recession_length,'n_start',n_start);
+
 end
 
 % add results to struct array
@@ -195,7 +201,9 @@
 results.Recession_a_Seasonality_error_str = Recession_a_Seasonality_error_str;
 results.AverageStorage = AverageStorage;
 results.AverageStorage_error_str = AverageStorage_error_str;
-results.RecessionParameters = RecessionParameters;
+results.RecessionParameters_a = RecessionParameters_a;
+results.RecessionParameters_b = RecessionParameters_b;
+results.RecessionParameters_T0 = RecessionParameters_T0;
 results.RecessionParameters_error_str = RecessionParameters_error_str;
 results.MRC_num_segments = MRC_num_segments;
 results.MRC_num_segments_error_str = MRC_num_segments_error_str;

TOSSH_code/signature_functions/sig_RecessionAnalysis.m

@@ -25,7 +25,9 @@
 %
 %   OUTPUT
 %   Recession_Parameters: matrix with parameters alpha, beta (=1 for 
-%       exponential fit) for each recession segment
+%       exponential fit) for each recession segment.
+%       Recession_Parameters(:,1): a, scaling parameter
+%       Recession_Parameters(:,2): b, parameter of non-linearity
 %   recession_month: approx. month of recession
 %   error_flag: 0 (no error), 1 (warning), 2 (error in data check), 3
 %       (error in signature calculation)
@@ -57,6 +59,37 @@
 %   This software is distributed under the GNU Public License Version 3.
 %   See <https://www.gnu.org/licenses/gpl-3.0.en.html> for details.
 
+%   NOTE: Post-processing of Recession Parameters
+%
+%   Post-processing of Recession parameters can be done as:
+%   RecessionParameters_a(i) = median((RecessionParametersTemp(:,1)),'omitnan');
+%   RecessionParameters_b(i) = median(RecessionParametersTemp(:,2),'omitnan');
+%   RecessionParametersT0Temp = 1./(RecessionParametersTemp(:,1).*median(Q_mat{i}(Q_mat{i}>0),'omitnan').^(RecessionParametersTemp(:,2)-1));
+%   ReasonableT0 = and(RecessionParametersTemp(:,2)>0.5,RecessionParametersTemp(:,2)<5);
+%   RecessionParameters_T0(i) = median(RecessionParametersT0Temp(ReasonableT0),'omitnan');
+
+%   RecessionParameters_T0 is characteristic timescale of recessions at median flow; It can be obtained 
+%   by fitting a line to the dQ/dt versus Q point cloud in log-log space for each individual recession, 
+%   with  Q scaled by median Q; T0 is the median value of −1/intercept (McMillan et al., 2021). 
+
+%   RecessionParameters_T0 can be derived from RecessionParameters_a and RecessionParameters_b as follows (McMillan et al, 2014): 
+%   Let Qhat = Q/Qmedian, then:
+%       dQhat/dt = - (1/T0) * Qhat^b
+%   Separating constant terms,
+%       1/Qmedian * dQ/dt = - (1/T0) * Q^b * (1/Qmedian)^b
+%       dQ/dt = - (1/T0) * Q^b * (1/Qmedian) ^ (b-1)
+%   As our estimate of a and b parameters are from:
+%       dQ/dt = - a * Q^b
+%   It leads to: 
+%       a = (1/T0) * (1/Qmedian) ^ (b-1)
+%       T0 = (1/a) * (1/Qmedian) ^ (b-1)
+%       T0 = 1 / (a * Qmedian ^ (b-1)) # This corresponds to the code for calculating RecessionParametersT0Temp
+%
+%   McMillan, H., Gueguen, M., Grimon, E., Woods, R., Clark, M., & Rupp, D. E. (2014). 
+%   Spatial variability of hydrological processes and model structure diagnostics in a 50 km2 catchment. 
+%   Hydrological Processes, 28(18), 4896-4913. https://doi.org/10.1002/hyp.9988
+
+
 % check input parameters
 if nargin < 2
     error('Not enough input arguments.')