1 /********************************************************** 2 * Version $Id: ihacres_eq.h 911 2011-02-14 16:38:15Z reklov_w $ 3 *********************************************************/ 4 /////////////////////////////////////////////////////////// 5 // ihacres_eq.h // 6 // // 7 // Copyright (C) 2006 by // 8 // Stefan Liersch // 9 //-------------------------------------------------------// 10 // e-mail: stefan.liersch@ufz.de // 11 // stefan.liersch@gmail.com // 12 // 2006-08-27 // 13 // modified: 2008-01-28 // 14 //-------------------------------------------------------// 15 // 16 // References: 17 // Jakeman, A.J. / Hornberger, G.M (1993). 18 // How Much Complexity Is Warranted in a 19 // Rainfall-Runoff Model? 20 // Water Resources Research, (29), NO. 8 (2637-2649) 21 // Kokkonen, T. S. et al. (2003). 22 // Predicting daily flows in ungauged catchments: 23 // model regionalization from catchment descriptors 24 // at the Coweeta Hydrologic Laboratory, North Carolina 25 // Hydrological Processes (17), 2219-2238 26 // Croke, B. F. W., W. S. Merritt, et al. (2004). 27 // A dynamic model for predicting hydrologic response 28 // to land cover changes in gauged and 29 // ungauged catchments. 30 // Journal Of Hydrology 291(1-2): 115-131. 31 // Croke, B. F. W., Andrews, F., Jakeman, A. J., Cuddy, S., Luddy, A.: 32 // (2005). Redesign of the IHACRES rainfall-runoff model 33 /////////////////////////////////////////////////////////// 34 35 //*******************************************************// 36 // NOTES // 37 //*******************************************************// 38 // This class provides functions of the rainfall-runoff 39 // model IHACRES. On the one hand the class can be used as 40 // a provider of functions using the standard constructor 41 // or the class can be instantiated and then performs the 42 // calculations automatically. 43 // 44 // Every function working with fields can be provided with 45 // arrays (pointer) or vector. You find for both solutions 46 // the appropriate function. 47 //*******************************************************// 48 49 //*******************************************************// 50 // ToDo // 51 //-------------------------------------------------------// 52 // - Difference eRain-streamflow_obs 53 // - Tau(q) and Tau(s) 54 // - Tau-values for different storages 55 // - if TMP_data_exist = false ??? 56 // - hourly data (not only daily) 57 // 58 // - SNOW TOOL 59 // - snow module on/off 60 // - snow module integration in CalcExcessRain (add bool variable) 61 // - implementation of an instanteneous store 62 // - implement more efficiency criterions 63 //*******************************************************// 64 65 //--------------------------------------------------------- 66 #ifndef HEADER_INCLUDED__ihacres_eq_H 67 #define HEADER_INCLUDED__ihacres_eq_H 68 //--------------------------------------------------------- 69 70 #include <saga_api/saga_api.h> // CSG_Table 71 72 #include <vector> // used for storing date string values in array 73 // and double arrays (streamflow, pcp, tmp ...) 74 #include <string> 75 76 using namespace std; 77 78 #include "snow_module.h" 79 #include "model_tools.h" 80 81 typedef std::vector<std::string> date_array; 82 typedef vector<double> vector_d; 83 //--------------------------------------------------------- 84 85 //--------------------------------------------------------- 86 // A CLASS TO DEAL WITH PARAMETERS OF THE 87 // LINEAR STORAGE TOOL 88 //--------------------------------------------------------- 89 class C_IHAC_LinearParms 90 { 91 public: C_IHAC_LinearParms()92 C_IHAC_LinearParms() { 93 _ZeroPointers(); 94 } C_IHAC_LinearParms(int size,int nStorages)95 C_IHAC_LinearParms(int size, int nStorages) { 96 _ZeroPointers(); 97 this->nStorages = nStorages; 98 if (nStorages == 1) 99 { 100 a = new double[size]; 101 b = new double[size]; 102 } 103 if (nStorages == 2) 104 { 105 aq = new double[size]; 106 as = new double[size]; 107 bq = new double[size]; 108 bs = new double[size]; 109 } 110 } ~C_IHAC_LinearParms(void)111 ~C_IHAC_LinearParms(void) { 112 if (nStorages == 1) 113 { 114 if (a) delete[] a; 115 if (b) delete[] b; 116 } 117 if (nStorages == 2) 118 { 119 if (aq) delete[] aq; 120 if (as) delete[] as; 121 if (bq) delete[] bq; 122 if (bs) delete[] bs; 123 } 124 } 125 126 // linear module parameters 127 int nStorages; 128 double* a; 129 double* b; 130 double* aq; 131 double* as; 132 double* bq; 133 double* bs; 134 135 private: _ZeroPointers()136 void _ZeroPointers() { 137 a = NULL; 138 b = NULL; 139 aq = NULL; 140 as = NULL; 141 bq = NULL; 142 bs = NULL; 143 } 144 }; 145 //--------------------------------------------------------- 146 147 //--------------------------------------------------------- 148 // A CLASS TO DEAL WITH PARAMETERS OF THE 149 // NON-LINEAR RAINFALL-LOSSES TOOL 150 //--------------------------------------------------------- 151 class C_IHAC_NonLinearParms 152 { 153 public: C_IHAC_NonLinearParms()154 C_IHAC_NonLinearParms() { 155 mp_tw = NULL; 156 mp_f = NULL; 157 mp_c = NULL; 158 mp_l = NULL; 159 mp_p = NULL; 160 mp_eR_flow_dif = NULL; 161 }; C_IHAC_NonLinearParms(int size)162 C_IHAC_NonLinearParms(int size) { 163 mp_tw = new double[size]; 164 mp_f = new double[size]; 165 mp_c = new double[size]; 166 mp_l = new double[size]; 167 mp_p = new double[size]; 168 mp_eR_flow_dif = new double[size]; 169 }; ~C_IHAC_NonLinearParms(void)170 ~C_IHAC_NonLinearParms(void) { 171 if (mp_tw) delete[] mp_tw; 172 if (mp_f) delete[] mp_f; 173 if (mp_c) delete[] mp_c; 174 if (mp_l) delete[] mp_l; 175 if (mp_p) delete[] mp_p; 176 if (mp_eR_flow_dif) delete[] mp_eR_flow_dif; 177 }; 178 // non-linear module parms 179 double* mp_tw; 180 double* mp_f; 181 double* mp_c; 182 double* mp_l; 183 double* mp_p; 184 double* mp_eR_flow_dif; 185 }; 186 //--------------------------------------------------------- 187 188 189 /////////////////////////////////////////////////////////////////////// 190 // 191 // CLASS Cihacres_eq 192 // 193 /////////////////////////////////////////////////////////////////////// 194 class Cihacres_eq 195 { 196 public: 197 198 /////////////////////////////////////////////////////////////////// 199 // 200 // CONSTRUCTORS 201 // 202 /////////////////////////////////////////////////////////////////// 203 204 // default 205 Cihacres_eq(); 206 207 // two storages 208 // using vector<double> as input 209 // if no temperature data are available 210 Cihacres_eq (date_array date, 211 vector_d streamflow_obs, 212 vector_d precipitation, 213 double Tw, double f, double c, 214 double l, double p, 215 double aq, double as, double bq, double bs); 216 // two storages 217 // using vector<double> as input 218 // if temperature data are available 219 Cihacres_eq (date_array date, 220 vector_d streamflow_obs, 221 vector_d precipitation, 222 vector_d temperature, 223 double Tw, double f, double c, 224 double l, double p, 225 double aq, double as, double bq, double bs, 226 double area, bool TMP_data_exist, 227 int IHAC_vers, 228 int storconf, 229 bool bSnowModule, 230 CSnowModule *SnowMod, 231 //double T_Rain, 232 //double T_Melt, 233 //double DD_FAC, 234 int delay); 235 236 // two storages 237 // using double arrays as input 238 // if no temperature data are available 239 Cihacres_eq (int size, // array size 240 date_array date, 241 double *streamflow_obs, 242 double *precipitation, 243 double Tw, double f, double c, 244 double aq, double as, double bq, double bs); 245 // two storages 246 // if temperature data are available 247 // using double arrays as input 248 Cihacres_eq (int size, // array size 249 date_array date, 250 double *streamflow_obs, 251 double *precipitation, 252 double *temperature, 253 double Tw, double f, double c, 254 double aq, double as, double bq, double bs); 255 256 // end constructors /////////////////////////////////////////////// 257 258 // destructor 259 ~Cihacres_eq(void); 260 261 /////////////////////////////////////////////////////////////////// 262 // 263 // PUBLIC MEMBER VARIABLES 264 // 265 /////////////////////////////////////////////////////////////////// 266 267 // Variables 268 int sizeAll; // incoming data array size (number of records in input data) 269 270 CSnowModule *m_pSnowMod; // class Snow Tool 271 272 double* m_pMeltRate; 273 274 /////////////////////////////////////////////////////////////////// 275 // 276 // PUBLIC MEMBER FUNCTIONS 277 // 278 /////////////////////////////////////////////////////////////////// 279 280 //-------------------------------------------------------- 281 282 void RunNonLinearModule (bool TMP_data_exist, bool bSnowModule, double T_Rain); 283 284 // simulate streamflow (single storage) 285 void SimStreamflowSingle (vector_d &excessRain, double initVal, 286 vector_d &streamflow_sim, int delay, 287 double a, double b); 288 289 void SimStreamflowSingle (double *excessRain, double initVal, 290 double *streamflow_sim, int delay, 291 double a, double b, int size); 292 293 // simulate streamflow with two parallel storages 294 void SimStreamflow2Parallel (vector_d &excessRain, 295 vector_d &streamflow_sim, 296 double initVal, // first observed streamflow value 297 double aq, double as, double bq, double bs, 298 double &vq, double &vs, 299 int IHAC_vers, 300 int delay); 301 302 void SimStreamflow2Parallel (double *excessRain, 303 double *streamflow_sim, 304 double initVal, // first observed streamflow value 305 double aq, double as, double bq, double bs, 306 double &vq, double &vs, 307 int IHAC_vers, int size, 308 int delay); 309 310 void SimStreamflow2Parallel (double *excessRain, 311 double *streamflow_sim, 312 double initVal, // first observed streamflow value 313 C_IHAC_LinearParms* linparms, int index, 314 double &vq, double &vs, int size, int delay); 315 316 double Calc_Parm_BS (double aq, double as, double bq); 317 318 // calculate time of decay for quick or slow component (aq or bq) 319 double Calc_TimeOfDecay (double a); 320 321 //-------------------------------------------------------- 322 323 // calculating the time constant, or inversely the rate at which 324 // the catchment wetness declines in the absence of rainfall 325 void CalcWetnessTimeConst (vector_d &temperature, 326 vector_d &Tw, 327 double TwConst, 328 double f); 329 330 void CalcWetnessTimeConst (double *temperature, 331 double *Tw, 332 double TwConst, 333 double f, 334 int size); 335 336 void CalcWetnessTimeConst (double* temperature, double* Tw, 337 C_IHAC_NonLinearParms* nonlinparms, int index, 338 int size); 339 340 // For ihacres_climate_scen 341 void CalcWetnessTimeConst_scen(double* temperature, double* Tw, 342 C_IHAC_NonLinearParms* nonlinparms, int index, 343 int size); 344 345 // modified version after Croke et al. (2005) 346 void CalcWetnessTimeConst_Redesign(vector_d &temperature, 347 vector_d &Tw, 348 double TwConst, 349 double f); 350 351 // modified version after Croke et al. (2005) 352 void CalcWetnessTimeConst_Redesign(double *temperature, 353 double *Tw, 354 double TwConst, 355 double f, 356 int size); 357 358 void CalcWetnessTimeConst_Redesign(double *temperature, double *Tw, 359 C_IHAC_NonLinearParms* nonlinparms, int index, 360 int size); 361 362 // calculating the catchment wetness index, 363 // or antecedent precipitation index 364 void CalcWetnessIndex (vector_d &Tw, 365 vector_d &precipitation, vector_d &temperature, 366 vector_d &WetnessIndex, double WI_init, 367 double c, bool bSnowModule, double T_Rain); 368 369 void CalcWetnessIndex (double *Tw, 370 double *precipitation, double *temperature, 371 double *WetnessIndex, double WI_init, 372 double c, bool bSnowModule, double T_Rain, 373 int size); 374 375 // modified version after Croke et al. (2005) 376 void CalcWetnessIndex_Redesign(vector_d &Tw, 377 vector_d &precipitation, 378 vector_d &WetnessIndex, 379 bool bSnowModule,double T_Rain); 380 381 // modified version after Croke et al. (2005) 382 void CalcWetnessIndex_Redesign(double *Tw, 383 double *precipitation, 384 double *WetnessIndex, double WI_init, 385 bool bSnowModule,double T_Rain, 386 int size); 387 388 // calculate the effecive or excess rainfall and 389 // calculate total effective rainfall over the period in [mm] 390 double CalcExcessRain (vector_d &precipitation, 391 vector_d &temperature, 392 vector_d &WetnessIndex, 393 vector_d &excessRain, double eR_init, 394 double &sum_eRainGTpcp, 395 bool bSnowModule, 396 CSnowModule* pSnowModule); 397 398 double CalcExcessRain (double *precipitation, 399 double *temperature, 400 double *WetnessIndex, 401 double *excessRain, double eR_init, 402 double &sum_eRainGTpcp, 403 int size, 404 bool bSnowModule, double T_Rain, double T_Melt, double* MeltRate); 405 406 // modified version after Croke et al. (2005) 407 double CalcExcessRain_Redesign (vector_d &precipitation, 408 vector_d &temperature, 409 vector_d &WetnessIndex, 410 vector_d &excessRain, double eR_init, 411 double &sum_eRainGTpcp, 412 double c, // mass balance parameter 413 double l, // soil moisture index threshold 414 double p, // power on soil moisture 415 bool bSnowModule, CSnowModule* m_pSnowMod); 416 417 // modified version after Croke et al. (2005) 418 double CalcExcessRain_Redesign (double *precipitation, 419 double *temperature, 420 double *WetnessIndex, 421 double *excessRain, double eR_init, 422 double &sum_eRainGTpcp, 423 int size, 424 double c, // mass balance parameter 425 double l, // soil moisture index threshold 426 double p, // power on soil moisture 427 bool bSnowModule, double T_Rain, double T_Melt, double* MeltRate); 428 //-------------------------------------------------------- 429 430 431 //-------------------------------------------------------- 432 void AssignFirstLastRec (CSG_Table &pTable, 433 int &first, int &last, 434 CSG_String date1,CSG_String date2, 435 int dateField); 436 437 int Assign_nStorages (int storconf); 438 439 double SumVector (vector_d &input); 440 double SumVector (double *input, int size); 441 442 double _Assign_NSE_temp (int obj_func, 443 double NSE, 444 double NSE_highflow, 445 double NSE_lowflow); 446 447 //-------------------------------------------------------- 448 // Get Functions 449 //-------------------------------------------------------- 450 double get_vq(); // fraction of quick flow 451 double get_vs(); // fraction of slow flow 452 double get_sum_streamflowMM_obs(int size); // sum of observed streamflow in [mm] 453 double get_sum_precipitation(int size); 454 double get_NSE(); 455 456 vector_d get_streamflow_sim(); 457 vector_d get_excessRain(); 458 vector_d get_WetnessIndex(); 459 vector_d get_Tw(); 460 461 private: 462 463 /////////////////////////////////////////////////////////////////// 464 // Private Member Variables 465 /////////////////////////////////////////////////////////////////// 466 467 // incoming arrays / vectors 468 date_array date; // Vector containing date values 469 vector_d streamflow_obs; // Vector containing observed streamflow data 470 vector_d precipitation; // Vector containing measured precipitation 471 vector_d temperature; // Vector containing measured temperature 472 vector_d streamflowMM_obs; // Streamflow time series in [mm] 473 474 // 475 double sum_eRainGTpcp; // sum of excess rainfall that is greater 476 // than precipitation in one simulation 477 478 // Vectors (arrays) 479 vector_d streamflow_sim; // Vector containing simulated streamflow data 480 vector_d excessRain; // Vector containing estimated excess or effective rainfall 481 vector_d WetnessIndex; // catchment wetness index 482 vector_d Tw; 483 484 // parameters 485 bool TMP_data_exist; // if temperature data are available (TMP_data_exist = true) 486 double RainRunoffCoef; // Rainfall-Runoff coefficient 487 double NSE; // Nash-Sutcliffe Efficiency 488 489 double sum_streamflow_obs; 490 double sum_streamflow_sim; 491 double sum_eRainMM; // sum of estimated effective rainfall [mm] 492 493 // non-linear IHACRES module parameters 494 double c; 495 double f; 496 double TwConst; // Tw is approximately the time constant, or inversely, 497 // the rate at which the catchment wetness declines 498 // in the absence of rainfall. 499 // Additional parameters for Croke et al. (2005) Redesign version 500 double l; // soil moisture index threshold 501 double p; // power on soil moisture 502 503 // linear IHACRES module parameters 504 double a; 505 double b; 506 507 double aq; 508 double as; 509 double bq; 510 double bs; 511 512 double vq; // fraction of quick flow 513 double vs; // fraction of slow flow 514 515 int delay; // The delay after the start of rainfall, 516 // before the discharge starts to rise. 517 518 double area; // area of the watershed in [km2] 519 520 int IHAC_version; // Different versions of IHACRES exist, corresponding 521 // to the version... 522 523 bool bSnowModule; // if true, snow module is active 524 525 /* 526 // snow module parameters 527 vector_d SnowStorage; 528 vector_d MeltRate; 529 double T_Rain; 530 double T_Melt; 531 double DD_FAC; 532 */ 533 /////////////////////////////////////////////////////////////////// 534 // Private Member Functions 535 /////////////////////////////////////////////////////////////////// 536 537 // initialize Vectors 538 void _InitVectorsStart (int size); 539 540 // Resize all vectors *.resize(0) 541 void _ZeroAllVectors (); 542 543 }; // end class ihacres_eq 544 545 //--------------------------------------------------------- 546 #endif /* #ifndef HEADER_INCLUDED__ihacres_eq_H */ 547