GR4J adjusting X2 and X3 for subbasins

[jimstagge]

Is it possible to vary the X2 and X3 parameters for the GR4J model across subbasins in Raven?

I see there are GR4J_X2 and GR4J_X3 options under :SoilParameterList, but because of the way Raven is mimics the GR4J model, it seems these options manually change X2 and X3 at the soil unit level (clay, loam, etc..), not at the lumped soil profile level.

It appears possible to adjust X1, via the SOIL_PROD at the profile level, and X4, by making dummy land use classes and then manually adjusting via LandUseParameterList .

I quickly modified the GR4J Salmon tutorial example:


# - Soil classes ---------------------------------------------
:SoilClasses
:Attributes
:Units
SOIL_PROD
SOIL_ROUT
SOIL_TEMP
SOIL_GW
:EndSoilClasses

# ----Soil Profiles--------------------------------------------
# name,#horizons,{soiltype,thickness}x{#horizons}
# GR4J_X1 is thickness of first layer (SOIL_PROD), here 0.529
:SoilProfiles
DEFAULT_P, 4, SOIL_PROD , 0.529, SOIL_ROUT , 0.300, SOIL_TEMP , 1.000, SOIL_GW , 1.000,
soil1, 4, SOIL_PROD , 0.700, SOIL_ROUT , 0.300, SOIL_TEMP , 1.000, SOIL_GW , 1.000,
soil2, 4, SOIL_PROD , 0.200, SOIL_ROUT , 0.300, SOIL_TEMP , 1.000, SOIL_GW , 1.000,
:EndSoilProfiles

# ----Vegetation Classes---------------------------------------
:VegetationClasses
:Attributes, MAX_HT, MAX_LAI, MAX_LEAF_COND
:Units, m, none, mm_per_s
VEG_ALL, 0.0, 0.0, 0.0
:EndVegetationClasses

# --Land Use Classes------------------------------------------
:LandUseClasses
:Attributes, IMPERM, FOREST_COV
:Units , frac, frac
LU_ALL, 0.0, 0.0
lu1, 0.0, 0.0
lu2, 0.0, 0.0

:EndLandUseClasses

# -Global Parameters -----------------------------------------
:GlobalParameter RAINSNOW_TEMP 0.0
:GlobalParameter RAINSNOW_DELTA 1.0
:GlobalParameter AIRSNOW_COEFF 0.053 # [1/d] = 1.0 - CEMANEIGE_X2 = 1.0 - x6
:GlobalParameter AVG_ANNUAL_SNOW 16.9 # [mm] = CEMANEIGE_X1 = x5
:GlobalParameter PRECIP_LAPSE 0.0004
:GlobalParameter ADIABATIC_LAPSE 0.0065
:GlobalParameter AVG_ANNUAL_RUNOFF 50

:SoilParameterList
:Parameters, POROSITY , GR4J_X3, GR4J_X2
:Units , none , mm, mm/d
[DEFAULT], 1.0 , 407.29, -3.396
soil1, 1.0 , 200, 3.396
soil2, 1.0 , 600, -3.396
:EndSoilParameterList

:LandUseParameterList
:Parameters, GR4J_X4, MELT_FACTOR
:Units , d, mm/d/C
[DEFAULT], 1.072, 7.73
lu1, 1, 7.73
lu2, 2, 7.73
:EndLandUseParameterList


Any help would be appreciated. For now I'm just runnning gr4j outside raven and dumping flows, but it would be nice to have everything running internally - thanks!

[rchlumsk]

Hi Jim

Apologies for the slow response. This is possible to do in Raven - I have provided an example below, further modifying the Salmon example. The main thing is to know which soil class each parameter is associated with, which is documented in the manual. In your example the X4 and X1 parameters can also be modified, included here as well. The example shows the X2 and X3 parameters both defined for the ROUTING_STORE and TEMP_STORE layers in percolation, and these are modified in the rvp file for subbasins 1 and 2. The example below is with 2 subbasins and 2 hrus.

Code: Select all

#########################################################################                                 
:FileType          rvp ASCII Raven 2.8.2                                                                             
#
# Emulation of GR4J simulation of Salmon River near Prince George                                                             
#------------------------------------------------------------------------

# - Soil classes ---------------------------------------------
:SoilClasses
  :Attributes
  :Units     
   SOIL_PROD1     
   SOIL_ROUT1
   SOIL_TEMP1
   SOIL_GW1
   SOIL_PROD2   
   SOIL_ROUT2
   SOIL_TEMP2
   SOIL_GW2   
:EndSoilClasses

# ----Soil Profiles--------------------------------------------
#     name,#horizons,{soiltype,thickness}x{#horizons}
#     GR4J_X1 is thickness of first layer (SOIL_PROD), here 0.529 and 0.500
:SoilProfiles
  DEFAULT_P1,      4, SOIL_PROD1   ,   0.529, SOIL_ROUT1  ,   0.300, SOIL_TEMP1  ,   1.000, SOIL_GW1  ,   1.000,
  DEFAULT_P2,      4, SOIL_PROD2   ,   0.500, SOIL_ROUT2  ,   0.300, SOIL_TEMP2  ,   1.000, SOIL_GW2  ,   1.000,
:EndSoilProfiles

# ----Vegetation Classes---------------------------------------
:VegetationClasses
   :Attributes,       MAX_HT,       MAX_LAI,      MAX_LEAF_COND
        :Units,            m,          none,           mm_per_s
       VEG_ALL,           0.0,          0.0,                0.0
:EndVegetationClasses

# --Land Use Classes------------------------------------------
:LandUseClasses
  :Attributes, IMPERM, FOREST_COV
  :Units     ,   frac,       frac
       LU_ALL,    0.0,        0.0
      LU1,         0.0,        0.0
      LU2,         0.0,        0.0
:EndLandUseClasses

# -Global Parameters  -----------------------------------------
:GlobalParameter RAINSNOW_TEMP       0.0
:GlobalParameter RAINSNOW_DELTA      1.0
:GlobalParameter AIRSNOW_COEFF     0.053 # [1/d] = 1.0 - CEMANEIGE_X2 = 1.0 - x6
:GlobalParameter AVG_ANNUAL_SNOW    16.9 # [mm]  =       CEMANEIGE_X1 =       x5
:GlobalParameter PRECIP_LAPSE     0.0004
:GlobalParameter ADIABATIC_LAPSE  0.0065

:SoilParameterList
 :Parameters, POROSITY ,  GR4J_X3, GR4J_X2
 :Units     ,     none ,       mm,    mm/d     
   [DEFAULT],      1.0 ,   407.29,  -3.396
   SOIL_TEMP1,     _DEFAULT,  401.29,  -2.396   # changing X2, X3 in subbasin 1
   SOIL_ROUT1,     _DEFAULT,  401.29,  -2.396   # changing X2, X3 in subbasin 1
   SOIL_TEMP2,     _DEFAULT,  427.29,  -4.396   # changing X2, X3 in subbasin 2
   SOIL_ROUT2,     _DEFAULT,  427.29,  -4.396   # changing X2, X3 in subbasin 2
:EndSoilParameterList

:LandUseParameterList
 :Parameters, GR4J_X4, MELT_FACTOR
 :Units     ,       d,      mm/d/C     
   [DEFAULT],   1.072,        7.73
   LU1,            1.00,     _DEFAULT   # changing X4 in subbasin 1
   LU2,            1.20,     _DEFAULT   # changing X4 in subbasin 2
:EndLandUseParameterList

:AvgAnnualRunoff 123.45  # placeholder average annual runoff, required with >1 subbasin


The rvh file I have is this:

Code: Select all

:SubBasins
  :Attributes,          NAME, DOWNSTREAM_ID,PROFILE,REACH_LENGTH,GAUGED
  :Units     ,          none,          none,   none,          km,  none
            1,   gr4j-salmon1,            -1,   NONE,       _AUTO,     1
         2,   gr4j-salmon2,            1,    NONE,       _AUTO,     1
:EndSubBasins
:HRUs
  :Attributes,  AREA, ELEVATION, LATITUDE, LONGITUDE, BASIN_ID,LAND_USE_CLASS, VEG_CLASS,SOIL_PROFILE, AQUIFER_PROFILE, TERRAIN_CLASS, SLOPE, ASPECT
  :Units     ,   km2,         m,      deg,       deg,     none,          none,      none,        none,            none,          none,   deg,    deg
            1,4250.6,     843.0,  54.4848, -123.3659,        1,        LU1,   VEG_ALL,   DEFAULT_P1,          [NONE],        [NONE],     0,      0
            2,4250.6,     843.0,  54.4848, -123.3659,        2,        LU2,   VEG_ALL,   DEFAULT_P2,          [NONE],        [NONE],     0,      0
:EndHRUs

One note here is that GR4J works reasonably well as a lumped model, and if you turn it into a distributed model, the number of parameters will increase greatly since the parameters are not tied to physical characteristics of HRUs and thus cannot be generalized across HRUs, as is typically done in distributed models. Depending on your application I would recommend that you either stick to a lumped GR4J model or consider a different model structure for distributed modelling, but the example shows how it could be done with GR4J.

Hope that helps!
Rob

[RMarinao]

Hi Rob

Thanks for the answer given to Jim, it has been helpful to me :)

I understand that the GR4J emulation is exact, as stated in your 2020 paper, but what can be said about the CemaNeige emulation?

At first impression I see that the parameters are not exactly the same, since CemaNeige_X1, originally corresponds to weighting coefficient for snow pack thermal state [-], and in Raven's rvp file something like AverageAnnualSnow [mm] appears; on the other hand, CemaNeige_X2, corresponds originally to degree-day melt coefficient [mm/°C/d], but in Raven there is an equivalence with 1 - AirSnowCoeffcient [1 / d], which has a different measurement unit.

I would appreciate your clarification.

By the way, many thanks to the Raven team for the great work done.

Rodrigo

[jrcraig]

Hi Rodrigo,

The Cema Neige routine was ported from an existing code, rather than from the literature, so there might be some minor differences, but I believe the details of the original algorithm are effectively emulated.

The potential melt should be calculated using the degree day approach, so the CemaNiege_X2 would correspond to the Raven parameter MELT_FACTOR. In Cema Neige the melt temperature DD_MELT_TEMP should be fixed at zero degrees. MELT_FACTOR corresponds to X2.

The snowmelt rate is given as (0.9*snow_cov+0.1)*MELT_FACTOR*(T-T_f), where T is the temperature, T_f is zero, and snow cover is calculated as the SWE/avg_annual_snow. The AverageAnnualSnow (in [mm]) corresponds to the G_seuil in Valery's original thesis, and it can be seen that Vmin is hardcoded at 0.1. This is only applied if the snow temperature is at zero degrees, where the snow temperature is modified using the AirSnowCoefficient [1/d], which is equivalent in purpose to the (1-C_TG) coefficient from Valery's thesis, but scaled such that it supports multiple time steps rather than simply daily time steps. These (AirSnowCoefficient and MELT_FACTOR) are the 2 free parameters referred to in that original source material. However, the AverageAnnualSnow may in practice be deemed a third free parameter.