GROUND.TXT Discussion of new GROUND subroutines 2005/03/09 Reasons for designing new GROUND subroutines -------------------------------------------- 1. The former ground subroutines did not handle snow very well. There are grid cells in southern Alaska that contain both land ice and ground. The land ice is melting significantly while the ground fraction is accumulating snow indefinitely. Less than 25% of the former scheme's coding was related to snow and ice, but more than 75% in new routines. 2. The former scheme used an explicit time scheme for calculating the surface fluxes, but the sensible heat flux occasionally had large oscillations (> 1000 W/m^2) from half hour to half hour. The new scheme use a complicated implicit time scheme (as is now used for sea ice and land ice), which considers the changes in surface properties and second ground layer properties when calculating the surface fluxes. 3. The former scheme produced too much evaporation (and not enough runoff) at low and mid latitudes. This was based on flow measurements at river mouths (90% are too low at low and mid latitudes), although basin wide precipitation seemed reasonable (sometimes too low or too high). There was a parameter in the former scheme that increased surface runoff, but it produced poorer results in my C-grid model. The former model's Amazon runoff was less than 40% of the observed value. 4. The water movement functions, although based on literature, were too sharp and were inappropriate for large GCM grid cells. The new scheme uses much simpler functions with weaker gradients. 5. To proceed with future improvements to ground hydrology (dynamic vegetation and carbon cycle modeling), I think it is best that I fully understand the capabilities and deficiencies of the scheme I start with. 6. The new ground hydrology scheme requires less computer time. The new scheme uses 4 layers with thicknesses in the ground: dZG(L) = .0625, .25, 1, and 4 m. The number of layers and thicknesses can be modified, but it does require some work. Soil fractions file that is used outside the GCM ------------------------------------------------ I have a program from Frank Abramopoulos (which he may have received from Robin Webb) that he used to define soil fractions for each grid box of his 6 layer model. Each grid box received a volumetric fraction of five different types: sand, silt, clay, peat and rock (which sum to 1). I modified this program for my 4 layer model and called it SOIL5X4.CRE . There are a few things I did not like about the unmodified results of this program: The global average rock fractions for my layers would have been: 0, .143, .247, .808 . Thus, I would have been left globally with only 1.8 m of soil which I think is too little. Also, there were large changes in the soil fractions from one layer to the next in the same column. To correct these problems, I replaced the fractions of each of the three deepest layers by averaging the fraction of a layer with the layer above it. Now the global non-rock thickness is almost 3 m. Based on the vegetation file, I replaced some soil with peat in the first layer. The reason is that vegetaion, litter, and shallow roots insulate the ground thermally, and peat has a much lower thermal conductivity than do other soil types. If there are objections to using peat, we could define an organic soil type, but we need to define all of the necessary parmeters. I am also considering reducing the large frations of peat in Canada and northern Europe. Soil properties file that is derived from the soil fractions file ----------------------------------------------------------------- The offline program, GRDAT5X4.CRE, creates the file GRDAT5X4.GISS which is read by the GCM in subroutine INPUT. The program requires the following quantities as a function of soil type: WFV0ofN = Minimum Volumetric Water Fraction WFVMofN = Maximum Volumetric Water Fraction HCDofN (10^6 J/C*m^2) = Dry Heat Capacity per unit Volume CTDofN (W/C*m) = Dry Thermal Conductivity CTWofN (W/C*m) = Thermal Conductivity with maximum Water CTIofN (W/C*m) = Thermal Conductivity with maximum Ice HDCofN (.001 kg/s*m) = Constant part of Hydraulic Diffusivity HDVofN (.001 kg/s*m) = Varying part of Hydraulic Diffusivity HDGofN (.001 kg/s*m^2)= Gravitational pull of water Diffusion RUNofN (1/year*m) = Underground source RUNoff removed Soil: Sand Silt Clay Peat Rock WFV0ofN .000 .018 .024 .000 .000 WFVMofN .396 .531 .576 .850 .000 HCDofN 2.0 2.0 2.0 2.5 2.0 CTDofN .9005 .3379 .2954 .0586 2.9000 CTWofN 3.3266 1.3988 1.2999 .5266 2.9000 CTIofN 4.6457 2.2455 2.1663 1.0138 2.9000 HDCofN .10 .10 .10 .10 .00 HDVofN .40 .20 .10 .80 .00 HDGofN .48 .12 .03 .24 .00 RUNofN .03 .01 .003 .015 .00 There are different sources for the parameters in this table. The output file, GRDAT5X4.GISS essentially contains, for each ground grid box, the quantities weighted by soil fractions. Most variables are defined within each layer. RUN is defined only for the three deepest layers. The HDx variables are defined between layer edges using 1/4 soil weighting of the upper layer and 3/4 for the lower layer. WFV0 comes from the former hydrology scheme's water fraction below which there is no water movement. WFVM and HCD come from published references. The only effect of WFV0 in the GCM is too slow down temperature changes because water and ice have relatively large heat capacities and because this water must be frozen or melted when the temperature is crossing 0 C. Layers can stay at 0 C for a long time. CTD, CTW and CTI come from replacing the former ground hydrology scheme's thermal conductivities with linear functions of water and ice. The thermal conductivity of a layer in the GCM is: CTD + CTW*WTR/WTRM + CTI*ICE/WTRM HDC, HDV and HDG are numbers I chose from general knowledge, they are not based on published observations. I have a set of graphs showing vertical water movement based on these coefficients and also a set using the former scheme's water movement which are considerably different. The hydraulic diffusivity used in the GCM is: HDC + HDV*WTRF1*WTRF2, where WTRF is the fraction (0 to 1) of liquid water in ground pores. I chose the variation of RUN among different soil types from general knowledge. The resultant coefficients are tuned by a global constant to improve river flow comparisons between the GCM and observations. The amount of water removed each time step from layers 2, 3 and 4 is: RUNU = dt*RUN*StDvZ*(WTR-WTR0)*(1-AIRF). StDvZ (m) is the standard deviation of .5 by .5 degree topography in a grid cell. Scaling of soil properties performed in subroutine INPUT -------------------------------------------------------- Subrountine INPUT reads in the file GRDAT5X4.GISS containing the soil weighted quantities listed above. These input quantites are scaled according to the following formulas: WTR0 = WFV0in * dZG(L) * RHOW WTRM = WFVMin * dZG(L) * RHOW HCD = HCDin * dZG(L) * 1000000 CTD = CTDin / dZG(L) dCTcW = CTWin / dZG(L) * WTRM dCTcI = CTIin / dZG(L) * WTRM HDC = HDCin * DTSURF / 1000 * .5*[dZG(L)+dZG(L+1)] HDV = HDVin * DTSURF / 1000 * .5*[dZG(L)+dZG(L+1)] HDG = HDGin * DTSURF / 1000 RUNX = RUNin * DTSURF * (WTRM-WTR0) / 365*SDAY byWTR = 1 / (WTRM-WTR0) DTSURF= 1800 s = time step used by subroutines GROUND and GRSNOW SDAY = 86400 s = number of seconds in a day WTR0 (kg/m^2) = minimum water mass in ground WTRM (kg/m^2) = maximum water mass in ground HCD (J/C*m^2) = dry heat capacity CTD (W/C*m^2) = dry themal conduction dCTdW (W/C*kg) = change of themal conduction by water dCTdI (W/C*kg) = change of themal conduction by ice HDC (kg/m^2) = constant part of hydraulic conduction HDV (kg/m^2) = varying part of hydraulic conduction HDG (kg/m^2) = gravitaional pull of water flux RUNX (kg/m^2) = underground source runoff coefficient byWTR (m^2/kg) = reciprocal of water maximum minus minimum Prognostic variables used in the new Ground Hydrology scheme ------------------------------------------------------------ RSN = horizontal fraction of ground that is snow covered MGR0 (kg/m^2) = mass of ice above layer 1, part of thermal layer 1 MGRn (kg/m^2) = mass of liquid water or ice in layer 1 without snow MSN (kg/m^2) = mass of snow in snow layer MGS0 (kg/m^2) = mass of ice above layer 1, part of thermal layer 1 MGSn (kg/m^2) = mass of liquid water or ice in layer 1 beneath snow HGRn (J/m^2) = heat content of thermal layer n without snow HSN (J/m^2) = heat content of snow layer HGSn (J/m^2) = heat content of thermal layer n beneath snow The ground surface fraction is subdivided horizontally into snow-free and snow-covered fractions specified by the prognostic variable RSN. Snow-free fraction has four thermal layers, HGRn, whose interfaces are specified by fixed thicknesses in the ground: .0625, .25, 1 and 4 m. Each thermal layer a has prognostic mass variable, MGRn, containing water and ice in the ground. The first thermal layer also contains ice above the ground, MGR0. MGR0 increases by compacted snow, snow melt that freezes, rain that freezes, and dew that freezes. The snow-covered fraction contains a single snow layer, MSN and HSN, plus ice and ground variables (MGS0, MGSn, HGSn) that reside beneath the snow and are comparable to the snow-free variables. Currently, the snow-free and snow-covered fractions receive the same radiation, but this will be corrected in the near future. How the modified subroutines are called --------------------------------------- The former ground hydrology scheme was usually called once each half hour during which precipitation and evaporation were mixed together. In the new scheme, an hour's worth of precipitation is applied to ground variables (PRECGR) immediately after atmospheric condensation. The half hour "surface loop" is now programmed in MAIN, inside which the surface fluxes are calculated and applied to the atmosphere (SURFCE) and then applied to the ground (GROUND, GRSNOW). PRECGR: apply precipitation to ground variables ----------------------------------------------- Inside PRECGR, precipitation is applied to the snow-free and snow- covered fractions separately. Twice the amount of rain that falls (dRNdSN=2) compresses snow into ice. Melt water from snow or ice is added to rain. Rain is partitioned into water that penetrates into layer 1 (dWTR) and surface runoff (RUNS). If the ground temperature is below freezing, then rain and melted snow water does not penetrate but is frozen and is added to the ice layer. If the ground temperature is at or above freeezing, then rain and melt water either penetrates or is added to source runoff (is added to the lakes). XINFIL = 1. dWTR = XINFIL * (RAIN+MELT) * (WTRM-MGR) / (WTRM-WTR0) = = XINFIL * (RAIN+MELT) * (air fraction of ground pores) RUNS = RAIN+MELT - dWTR Over the snow-free fraction when no rain is frozen at the surface, the rain that come from large scale condensation uses the formulas above, while the rain that penetrates and comes from moist convection uses: dWTR = MCNFIL * RAIN(MC) * (WTRM-MGR) / (WTRM-WTR0) MCNFIL = 1. Rind has been suggested that the infiltration factors (XINFIL and MCNFIL) should be dependent upon ground crusting in order to simulate crusting. Hillel has stated that infiltration is not a linear function of the rain rate; all rain infiltrates for slow rain, but above an infiltrability threshold, excessive rain runs off. These ideas should be implemented. PRECGR is the only subroutine that is allowed to change RSN. It is also the only subroutine that can increase the snow thickness. The first snowfall that does not melt on a grid cell without snow, divides the grid cell into a snow-free fraction and a snow-covered fraction with a snow thickness of 25.5 kg/m^2 (SNONEW). Additional snowfall increases the horizontal snow fraction at a smaller rate. The grid cell will become totally snow covered when the thickness reaches 2*SNONEW. If the snow thickness decreases (due to melting or evaporation) to below 24 kg/m^2 (SNOMIN), then the snow cover contracts horizontally so that the reduced snow cover is again at a thickness of SNONEW. SURFCE: calculate surface fluxes and apply them to atmosphere ------------------------------------------------------------- SURFCE is made simpler by eliminating the "surface loop" and by using separate coding for each surface type. It is now longer with duplicate code. The snow-free ground section now uses the more complicated implicit time scheme that will eliminate oscillations in the surface fluxes. The transition through 0 C is made more smoothly. Potential evaporation over ground is still modelled by a drag law formula. The ratio of actual evaporation to potential evaporation (WET), and which layers to draw the water from, depends upon several factors. Normally, WET is the ratio of current available water divided by maximum available water [(MGR-WTR0)/(WTRM-WTR0)] of layer 1. But over vegetation, when the sun is up and the ground temperature is above 5 C, a horizontal fraction of evaporation can be drawn from plants: FROOT = FVEG*ROOTT*ROOTS, where: FVEG = the vegetated fraction of the grid cell, ROOTT = 1 - (TCG1-23)^2/18^2, ROOTS = SRH/400 [SRH (W/m^2) = absorbed solar radiation at surface]. Thus the used water ratio (actual evaporation divided by potential) is: WETused = WET*(1-FROOT) + WETR*FROOT . WETR is a water ratio using abundant water in the first three layers. Later, WETR could be based on vegetation type. Current simulations with this new ground hydrology are quite promising. The most egregious current error is that most ground layers are exceedingly dry. Organic litter should be used to reduce evaporation from layer 1. Under snow-free ground, the diffusion of heat between ground layers 1 and 2 uses the harmonic mean of the thermal conductivity of each layer: RHOI (kg/m^3) = 910 = density of ice ALAMI (W/C*m) = 2.1762 = thermal conductivity of ice CTG (W/C*m^2) = CTD + WTR*dCTdW + ICE*dCTdI = conduction through layer dF1dTG(W/C*m^2) = 2 / (MGR0/RHOI*ALAMI + 1/CTG1 + 1/CTG2) F1 (W/m^2) = (TCG1 - TCG2) * dF1dTG = explicit heat flux Under snow-covered ground, the diffusion of heat between the snow layer and ground layers 1 uses the harmonic mean of the thermal conductivity of each layer: RHOS (kg/m^3) = 100 = density of snow ALAMI (W/C*m) = .35 = thermal conductivity of snow CTG1 (W/C*m^2) = CTD + WTR*dCTdW + ICE*dCTdI = conduction through layr1 dF1dTG(W/C*m^2) = 2 / (MSN/RHOS*ALAMS + MGS0/RHOI*ALAMI + 1/CTG1) F1 (W/m^2) = (TSNOW - TCG1) * dF1dTG = explicit heat flux The actual heat flux, F1, is derived implicitly. GROUND: surface fluxes are applied to snow-free ground variables ---------------------------------------------------------------- GROUND performs diffusion of heat, checks for melting of the ice layer (which will cause some surface runoff), applies dew or evaporation, diffuses liquid water, and removes liquid water from the deeper layers into lakes. Diffusion of heat between two layers uses the harmonic mean of thermal conductivity of each layer: CTG (W/C*m^2) = CTD + WTR*dCTdW + ICE*dCTdI = conduction through layer dF2dTG(W/C*m^2) = 2 / (1/CTG2 + 1/CTG3) F2 (W/m^2) = (TCG2 - TCG3) * dF2dTG = explicit heat flux The diffusive water flux is: WFLUX (kg/m^2) = (HDC + HDV*WTRF1*WTRF2) * (WTRF1*AIRF2 - WTRF2*AIRF1) + HDG*WTRF1*AIRF2 where: HDC (kg/m^2) = constant part of hydraulic diffusion HDV (kg/m^2) = varying part of hydraulic diffusion HDG (kg/m^2) = gravitational pull of water diffusion WTRF = liquid water fraction in pores, 0 to 1 AIRF = air fraction in pores, 0 to 1 If there is no ice, then AIRF = 1 - WTRF , and: WFLUX = (HDC + HDV*WTRF1*WTRF2) * (WTRF1 - WTRF2) + HDG*WTRF1*AIRF2 David Rind suggests pulling water out of rivers and into the ground in dry areas. GRSNOW: surface fluxes are applied to snow-covered ground variables ------------------------------------------------------------------- GRSNOW performs diffusion of heat, checks for melting of snow or ice (which will cause some surface runoff), applies dew or evaporation, diffuses liquid water, and removes liquid water from the deeper layers into lakes. Frozen dew increases the ice layer, not the snow layer. Snow that is melted from the top (but not from below) compresses additional snow (MELT*dSNdML = MELT*2) into ice. Currently, MELT water penetrates or runs off from the ground beneath the snow. Those variables could become nearly saturated which would cause most additional water to runoff. Some of this melt water should have a chance to penetrate into the snow-free ground. If a grid cell has minimum snow amount, is totally snow-covered, and is going to melt entirely, then RSN of the melt water should be applied to the snow- covered variables and 1-RSN should be applied to the snow-free variables as the cell is melting. References: DeVries D A, 1966. Thermal properties of soils. In "Physics of Plant Environment", Van Wijk W R ed., North-Holland Co., Amsterdam, 210-235. Zobler L, 1986. A world soil file for global climate modeling. NASA Technical Memorandum 87802. Abramopoulos F, Rosenzweig C, and Choudhury B, 1988. Improved ground hydrology calculations for global climate models (GCMs): soil water movement and evapotranspiration. Journal of Climate, 1(9), 921-941.