From 5318a3183c0f1ca8e222e10c2d0775eaaea3877f Mon Sep 17 00:00:00 2001 From: RY4GIT <400mhs2@gmail.com> Date: Fri, 14 Mar 2025 13:21:13 -0700 Subject: [PATCH] Streamlined RecessionAnalysis signature output and added more description --- TOSSH_code/calculation_functions/calc_All.m | 19 +++++- .../calc_McMillan_Groundwater.m | 58 +++++++++++-------- .../sig_RecessionAnalysis.m | 35 ++++++++++- 3 files changed, 83 insertions(+), 29 deletions(-) diff --git a/TOSSH_code/calculation_functions/calc_All.m b/TOSSH_code/calculation_functions/calc_All.m index 59ae735..6ca8d66 100644 --- a/TOSSH_code/calculation_functions/calc_All.m +++ b/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; diff --git a/TOSSH_code/calculation_functions/calc_McMillan_Groundwater.m b/TOSSH_code/calculation_functions/calc_McMillan_Groundwater.m index 783a7a4..42cfabc 100644 --- a/TOSSH_code/calculation_functions/calc_McMillan_Groundwater.m +++ b/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,52 +98,54 @@ % 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}, ... @@ -151,7 +153,7 @@ [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); @@ -159,12 +161,16 @@ 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; diff --git a/TOSSH_code/signature_functions/sig_RecessionAnalysis.m b/TOSSH_code/signature_functions/sig_RecessionAnalysis.m index 625fefe..4218eb1 100644 --- a/TOSSH_code/signature_functions/sig_RecessionAnalysis.m +++ b/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 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.')