- Generated on Mon Feb 3 2020 19:58:04 for CLASSIC by
1.8.13
|
CLASSIC
Canadian Land Surface Scheme including Biogeochemical Cycles
|
Calculates overland flow; steps ahead pond and soil layer temperatures, and checks for freezing of the ponded water and freezing or thawing of liquid or frozen water in the soil layers. Adjusts ponded water temperature, soil layer temperatures and water stores accordingly. More...
Functions/Subroutines | |
| subroutine | waterupdates (TBAR, THLIQ, THICE, HCP, TPOND, ZPOND, TSNOW, ZSNOW, ALBSNO, RHOSNO, HCPSNO, TBASE, OVRFLW, TOVRFL, RUNOFF, TRUNOF, HMFG, HTC, HTCS, WTRS, WTRG, FI, TBARW, GZERO, G12, G23, GGEO, TA, WSNOW, TCTOP, TCBOT, GFLUX, ZPLIM, THPOR, THLMIN, HCPS, DELZW, DELZZ, DELZ, ISAND, IWF, IG, ILG, IL1, IL2, JL, N) |
Calculates overland flow; steps ahead pond and soil layer temperatures, and checks for freezing of the ponded water and freezing or thawing of liquid or frozen water in the soil layers. Adjusts ponded water temperature, soil layer temperatures and water stores accordingly.
| subroutine waterupdates | ( | real, dimension (ilg,ig), intent(inout) | TBAR, |
| real, dimension (ilg,ig), intent(inout) | THLIQ, | ||
| real, dimension (ilg,ig), intent(inout) | THICE, | ||
| real, dimension (ilg,ig), intent(inout) | HCP, | ||
| real, dimension (ilg), intent(inout) | TPOND, | ||
| real, dimension (ilg), intent(inout) | ZPOND, | ||
| real, dimension (ilg), intent(inout) | TSNOW, | ||
| real, dimension (ilg), intent(inout) | ZSNOW, | ||
| real, dimension(ilg), intent(out) | ALBSNO, | ||
| real, dimension(ilg), intent(inout) | RHOSNO, | ||
| real, dimension(ilg), intent(inout) | HCPSNO, | ||
| real, dimension (ilg), intent(inout) | TBASE, | ||
| real, dimension(ilg), intent(inout) | OVRFLW, | ||
| real, dimension(ilg), intent(inout) | TOVRFL, | ||
| real, dimension(ilg), intent(inout) | RUNOFF, | ||
| real, dimension(ilg), intent(inout) | TRUNOF, | ||
| real, dimension (ilg,ig), intent(inout) | HMFG, | ||
| real, dimension (ilg,ig), intent(inout) | HTC, | ||
| real, dimension (ilg), intent(inout) | HTCS, | ||
| real, dimension (ilg), intent(inout) | WTRS, | ||
| real, dimension (ilg), intent(inout) | WTRG, | ||
| real, dimension (ilg), intent(in) | FI, | ||
| real, dimension (ilg,ig), intent(in) | TBARW, | ||
| real, dimension (ilg), intent(inout) | GZERO, | ||
| real, dimension (ilg), intent(in) | G12, | ||
| real, dimension (ilg), intent(in) | G23, | ||
| real, dimension (ilg), intent(in) | GGEO, | ||
| real, dimension (ilg), intent(in) | TA, | ||
| real, dimension (ilg), intent(in) | WSNOW, | ||
| real, dimension (ilg,ig), intent(in) | TCTOP, | ||
| real, dimension (ilg,ig), intent(in) | TCBOT, | ||
| real, dimension (ilg,ig), intent(inout) | GFLUX, | ||
| real, dimension (ilg), intent(in) | ZPLIM, | ||
| real, dimension (ilg,ig), intent(in) | THPOR, | ||
| real, dimension(ilg,ig), intent(in) | THLMIN, | ||
| real, dimension (ilg,ig), intent(in) | HCPS, | ||
| real, dimension (ilg,ig), intent(in) | DELZW, | ||
| real, dimension (ilg,ig), intent(in) | DELZZ, | ||
| real, dimension (ig), intent(in) | DELZ, | ||
| integer, dimension (ilg,ig), intent(in) | ISAND, | ||
| integer, intent(in) | IWF, | ||
| integer, intent(in) | IG, | ||
| integer, intent(in) | ILG, | ||
| integer, intent(in) | IL1, | ||
| integer, intent(in) | IL2, | ||
| integer, intent(in) | JL, | ||
| integer, intent(in) | N | ||
| ) |
| [in,out] | tbar | Temperature of soil layer \([C] (T_g)\) |
| [in,out] | thliq | Volumetric liquid water content of soil layer \((\theta_l) [m^3 m^{-3}]\) |
| [in,out] | thice | Volumetric frozen water content of soil layer \((\theta_i) [m^3 m^{-3}]\) |
| [in,out] | hcp | Heat capacity of soil layer \([J m^{-3} K^{-1}] (C_g)\) |
| [in,out] | hmfg | Energy associated with freezing or thawing of water in soil layer \([W m^{-2}]\) |
| [in,out] | htc | Internal energy change of soil layer due to conduction and/or change in mass \([W m^{-2}] (I_g)\) |
| [in,out] | tpond | Temperature of ponded water \([C] (T_p)\) |
| [in,out] | zpond | Depth of ponded water \([m] (z_p)\) |
| [in,out] | tsnow | Temperature of the snow pack \([C] (T_s)\) |
| [in,out] | zsnow | Depth of snow pack \([m] (z_g)\) |
| [out] | albsno | Albedo of snow [ ] |
| [in,out] | rhosno | Density of snow pack \((\rho_s) [kg m^{-3}]\) |
| [in,out] | hcpsno | Heat capacity of snow pack \([J m^{-3} K^{-1}] (C_s)\) |
| [in,out] | tbase | Temperature of bedrock in third soil layer,if only three layers are being modelled [K] |
| [in,out] | ovrflw | Overland flow from top of soil column [m] |
| [in,out] | tovrfl | Temperature of overland flow from top of soil column [K] |
| [in,out] | runoff | Total runoff from soil column [m] |
| [in,out] | trunof | Temperature of total runoff from soil column [K] |
| [in,out] | htcs | Internal energy change of snow pack due to conduction and/or change in mass \([W m^{-2}] (I_s)\) |
| [in,out] | wtrs | Water transferred into or out of the snow pack \([kg m^{-2} s^{-1}]\) |
| [in,out] | wtrg | Water transferred into or out of the soil \([kg m^{-2} s^{-1}]\) |
| [in] | fi | Fractional coverage of subarea in question on modelled area \([ ] (X_i)\) |
| [in] | tbarw | Temperature of water in soil layer [C] |
| [in,out] | gzero | Heat flow into soil surface \([W m^{-2}]\) |
| [in] | g12 | Heat flow between first and second soil layers \([W m^{-2}]\) |
| [in] | g23 | Heat flow between second and third soil layers \([W m^{-2}]\) |
| [in] | ggeo | Geothermal heat flux at bottom of soil profile \([W m^{-2}]\) |
| [in] | ta | Air temperature [K] |
| [in] | wsnow | Liquid water content of snow pack \([kg m^{-2}] (w_s)\) |
| [in] | tctop | Thermal conductivity of soil at top of soil layer \([W m^{-1} K^{-1}]\) |
| [in] | tcbot | Thermal conductivity of soil at bottom of soil layer \([W m^{-1} K^{-1}]\) |
| [in] | zplim | Limiting depth of ponded water [m] |
| [in] | thpor | Pore volume in soil layer \((\theta_p) [m^3 m^{-3}]\) |
| [in] | thlmin | Residual soil liquid water content remaining after freezing or evaporation \((\theta_r) [m^3 m^{-3}]\) |
| [in] | hcps | Heat capacity of soil material \([J m^{-3} K^{-1}] (C_m)\) |
| [in] | delzw | Permeable thickness of soil layer \((\Delta z_{g,w}) [m]\) |
| [in] | delzz | Soil layer thicknesses to bottom of permeable depth for standard three-layer configuration, or to bottom of thermal depth for multiple layers \((\Delta z_{g,z}) [m]\) |
| [in] | delz | Overall thickness of soil layer [m] |
| [in] | isand | Sand content flag |
| [in,out] | gflux | Heat flow between soil layers \([W m^{-2}]\) |
If a full-scale hydrological modelling application is not being run, that is, if only vertical fluxes of energy and moisture are being modelled, the flag IWF is pre-set to zero. In this case, overland flow of water is treated using a simple approach: if the ponded depth of water on the soil surface ZPOND exceeds a limiting value ZPLIM (which varies by land surface type), the excess is assigned to overland flow. These calculations are carried out in loop 100, for all modelled areas except ice sheets (which are handled in subroutine iceSheetBalance). The overall runoff from the modelled area in question, RUNOFF, is incremented by the excess ponded water, and the overland flow for the whole grid cell OVRFLW is incremented by the product of the excess ponded water and the fractional area of the grid cell. The temperature of the overall runoff from the modelled area TRUNOF, and the temperature of the overland flow for the grid cell TOVRFL, are calculated as weighted averages over their previous values and the ponded water temperature TPOND. Finally, ZPOND is set to ZPLIM.
In loop 200, the calculation of the change in internal energy HTC within the soil layers due to movement of soil water (addressed in the preceding subroutines waterFlowInfiltrate and waterFlowNonInfiltrate) is completed. (The initial step was performed at the end of subroutine soilWaterPhaseChg.) The volumetric heat capacity of the soil layers HCP is recalculated as a weighted average over the volumetric contents of liquid water, frozen water and soil particles. The temperature TBAR of each soil layer is recalculated as a weighted average of the liquid water temperature TBARW resulting from soil water movement, and the previous soil layer temperature. As noted in the documentation for subroutine soilWaterPhaseChg., TBARW and HCP apply only over the permeable depth of the soil layer, DELZW, whereas TBAR applies over the whole depth, DELZZ, and is associated with HCPSND, the volumetric heat capacity assigned to rock, in the interval DELZZ- DELZW. (For the default three-layer version of CLASS, DELZZ=DELZW in the third soil layer; see the documentation of loop 550 below.)
In loop 300 the calculation of the change of internal energy within the first soil layer due to changes in the surface ponded water is completed (for diagnostic purposes, the first level of HTC includes the internal energy of both the ponded water and the first soil layer). (The initial step of this calculation was performed at the end of subroutine pondedWaterFreeze.) The flow of heat between the bottom of the ponded water and the top of the first soil layer, GP1, is calculated by assuming a linear variation with depth of the heat flux between the top of the pond, GZERO, and between the first and second soil layers, G12. Since the heat flux G(z) varies directly with the temperature gradient dT(z)/dz, it can be seen that this approach is consistent with the assumption made in the energyBudgetDriver subroutines that T(z) is a quadratic function of depth (and therefore that its first derivative is a linear function of depth). The temperature of the ponded water is stepped ahead using GZERO and GP1, and GZERO is reset to GP1 for use in the later soil temperature calculations.
In the 400 loop, a check is performed to ascertain whether the ponded water temperature has fallen below 0 C as a result of the above temperature change. If so, calculations are carried out analogous to those done in pondedWaterFreeze. A recalculation of the available energy in the snow pack, HTCS, is initiated. The available energy sink HADD is calculated from TPOND and ZPOND, and TPOND is set to 0 C. The amount of energy required to freeze the whole depth of ponded water is calculated as HCONV. An update of the first level value of HTC, reflecting the freezing of ponded water, is initiated. If HADD < HCONV, the available energy sink is sufficient to freeze only part of the ponded water. This depth ZFREZ is calculated from HADD, and subtracted from ZPOND. ZFREZ is then converted from a depth of water to a depth of ice, and is reassigned as part of the snow pack. The snow temperature, density, heat capacity and depth are recalculated as weighted averages of the pre-existing snow properties and those of the frozen water. If there was no snow on the ground to begin with, the snow albedo is initialized to its minimum background value of 0.50.
If HADD > HCONV, the available energy sink is sufficient to freeze the whole depth of ponded water and then to decrease the frozen water temperature to below 0 C. ZFREZ is evaluated as ZPOND converted to a depth of ice, and the remaining energy sink HADD-HCONV is used to calculated a test temperature TTEST of the newly formed ice. In order to avoid unphysical overshoots of the frozen water temperature, a limiting temperature TLIM is determined as the minimum of the snow temperature and the first level soil temperature if there is pre-existing snow, and as the minimum of the air temperature and the first level soil temperature if not. If TTEST < TLIM, the excess energy sink required to decrease the frozen water temperature from TLIM to TTEST is calculated and added to GZERO, and the new frozen water temperature is assigned as TLIM; otherwise it is assigned as TTEST. The energy sink used to cool the frozen water from 0 C to TLIM or TTEST is decremented from the first level of HTC, and the snow temperature is recalculated as the weighted average of the pre-existing snow temperature and the frozen water temperature. If there was no pre-existing snow, the snow albedo is initialized to the background value of 0.50. The density, heat capacity and depth of the snow are recalculated as above. Finally, the recalculations of the internal energy diagnostic variables HTC and HTCS are completed; and ZFREZ is used to update the diagnostic variable HMFG describing the energy associated with phase changes of water in soil layers, and the diagnostic variables WTRS and WTRG describing transfers of water into or out of the snow and soil respectively.
In the 500, 550 and 600 loops, the heat fluxes between soil layers and the updated temperatures for the soil layers are calculated according to the discretization approach that has been selected for the model run. It is assumed that the first two layer thicknesses are always set to 0.10 and 0.25 m respectively, but two possible options exist for the treatment of deeper layers. In the operational, "three-layer" configuration, a single third soil layer is modelled of thickness 3.75 m, with a user-specified permeable thickness DELZW. TBAR of the third layer is taken to apply to this permeable thickness, whereas a separate bedrock temperature, TBASE, is carried for the impermeable thickness, represented by DELZ-DELZW. This strategy is adopted so that temperature variations caused by melting and freezing of water are confined to the permeable thickness. In the "multi-layer" configuration, there can be any number of deeper layers, of thicknesses specified by the user (with the proviso that the thicknesses not be small enough to lead to numerical instability in the forward explicit time-stepping scheme). In this configuration it is deemed unnecessary to carry a separate calculation of TBASE.
In the energyBudgetDriver subroutines, the heat fluxes at the top of the first soil layer, and between the first and second and the second and third soil layers, were determined. In loop 500, the fluxes between the remaining soil layers are calculated for the multi-layer configuration, and assigned to the heat flux vector GFLUX. Since the temperature variation is increasingly damped with depth, a simple linearization of the temperature profile is used. The expression for the ground heat flux at depth z, G(z), which depends on the thermal conductivity \(\lambda (z)\) and the temperature gradient, is given as: \(G(z) = \lambda (z) dT(z)/dz\) The linearized form is written as: \(G_j = (\lambda_{j-1} + \lambda_j) (T_{j-1} - T_j) / (\Delta z_{j-1} - \Delta z_j)\) where G j is the heat flux at the top of layer j, and \(\lambda_j\) , \(T_j\) and \(\Delta z_j\) refer to the thermal conductivity, temperature and heat capacity of the layer.
In the 550 loop, the first and second soil layer temperatures are stepped ahead using GZERO, G12 and G23. (Recall that the calculated heat capacity HCP applies to the permeable thickness DELZW, and the heat capacity of rock HCPSND applies to the rest of the layer thickness, DELZZ-DELZW.) If the three- layer configuration is being used, then a check is carried out to ascertain whether the permeable thickness DELZZ of the third layer is greater than zero. If so, and if DELZZ is not effectively equal to DELZ, the heat flow G3B at the interface between the permeable thickness and the underlying bedrock is calculated using the linearized heat flow equation given above; TBAR is updated using the difference between G23 and G3B, and TBASE is updated using the difference between G3B and GGEO, the geothermal heat flux at the bottom of the soil profile. If DELZZ is equal to DELZ, the whole layer is permeable and TBAR is updated using the difference between G23 and GGEO. If DELZZ is zero, the whole layer is impermeable, and TBASE is updated using the difference between G23 and GGEO. HTC for the third layer is updated with the GGEO value. Finally, if the multi-layer configuration is being used, TBAR of the third layer is updated using G23, the flux at the top of the layer. For diagnostic purposes, the first three levels of the GFLUX vector are assigned as GZERO, G12 and G23 respectively.
In the 600 loop, the remaining soil layer temperature calculations are done for layers 3 and deeper, for the multi-layer configuration. At the beginning and end of the loop an updated calculation of HTC for the layer in question is initiated and completed respectively. For the third layer, TBAR is updated using the heat flux at the bottom of the layer; for the last layer, TBAR is updated using the difference between GFLUX at the top of the layer and GGEO at the bottom. For the intermediate layers, TBAR is calculated using the difference between the GFLUX values at the top and bottom of the layer.
In the 700 loop, checks are carried out to determine whether, as a result of the forward stepping of the temperature, TBAR has fallen below 0 C while liquid water still exists in the layer, or risen above 0 C while frozen water still exists in the layer. If either occurs, adjustments to the water content are performed analogous to those done in subroutine soilWaterPhaseChg. Again, at the beginning and end of the loop an updated calculation of HTC for each layer in turn is initiated and completed respectively.
If the soil layer temperature is less than 0 C and the volumetric liquid water content THLIQ of the layer is greater than the residual water content THLMIN, the water content THFREZ that can be frozen by the available energy sink is calculated from TBAR and the weighted average of HCP and HCPSND. If THFREZ \(\leq\) THLIQ - THLMIN, all of the available energy sink is used to freeze part of the liquid water content in the permeable part of the soil layer. The amount of energy involved is subtracted from HTC and added to HMFG. THFREZ is subtracted from THLIQ, converted to an ice volume and added to THICE. HCP is recalculated, and the layer temperature is set to 0 C. Otherwise, all of the liquid water content of the layer above THLMIN is converted to frozen water, and HMFG and HTC are recalculated to reflect this. HCP is recomputed, and the remaining energy sink HADD is applied to decreasing the temperature of the soil layer.
If the soil layer temperature is greater than 0 C and the volumetric ice content THICE of the layer is greater than zero, the ice content THMELT that can be melted by the available energy is calculated from TBAR and the weighted average of HCP and HCPSND. If THMELT \(\leq\) THICE, all of the available energy is used to melt part of the frozen water content of the permeable part of the layer. The amount of energy involved is subtracted from HTC and added to HMFG. THMELT is subtracted from THICE, converted to a liquid water volume and added to THLIQ; HCP is recalculated and the layer temperature is set to 0 C. Otherwise, all of the frozen water content of the layer is converted to liquid water, and HMFG and HTC are recalculated to reflect this. HCP is recomputed, and the remaining energy HADD is applied to increasing the temperature of the soil layer.
1.8.13