Changeset 21398

Timestamp:

11/18/16 07:54:23 (8 years ago)

Author:

bdef

Message:

NEW:changing transmitivity and storing coeficient for confined-unconfined ability

Location:

issm/trunk-jpl/src/c/analyses

Files:

• HydrologyDCInefficientAnalysis.cpp (modified) (10 diffs)
• HydrologyDCInefficientAnalysis.h (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

@@ -190 +190 @@
         IssmDouble  D_scalar,Jdet,dt;
         IssmDouble  sediment_transmitivity;
-        IssmDouble  transfer;
+        IssmDouble  prestep_head,base_elev;
+        IssmDouble  transfer,sediment_storing;
         IssmDouble *xyz_list  = NULL;
 
@@ -215 +216 @@
         Input* base_input     = basalelement->GetInput(BaseEnum);
 
-        IssmDouble sediment_storing = SedimentStoring(basalelement);
         /*Transfer related Inputs*/
         if(isefficientlayer){
@@ -226 +226 @@
                 gauss          -> GaussPoint(ig);
                 basalelement   -> JacobianDeterminant(&Jdet,xyz_list,gauss);
-                SedTrans_input -> GetInputValue(&sediment_transmitivity,gauss);
+
+                sediment_transmitivity = SedimentTransmitivity(basalelement,gauss,sed_head_input,base_input,SedTrans_input);
+                sediment_storing       = SedimentStoring(basalelement,gauss,sed_head_input,base_input);
+
                 /*Diffusivity*/
                 D_scalar=sediment_transmitivity*gauss->weight*Jdet;
@@ -292 +295 @@
         /*Intermediaries */
         bool       active_element,isefficientlayer;
-        IssmDouble dt,scalar;
+        IssmDouble dt,scalar,sediment_storing;
         IssmDouble water_head;
         IssmDouble water_load,transfer;
@@ -313 +316 @@
         basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
 
-        Input* sed_head_input                     = basalelement->GetInput(SedimentHeadEnum);
-        Input* epl_head_input                     = basalelement->GetInput(EplHeadEnum);
-        Input* thick_input                                      = basalelement->GetInput(ThicknessEnum);
-        Input* base_input                                         = basalelement->GetInput(BaseEnum);
-        Input* water_input                                      = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input);
+        Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum);
+        Input* epl_head_input = basalelement->GetInput(EplHeadEnum);
+        Input* thick_input        = basalelement->GetInput(ThicknessEnum);
+        Input* base_input                 = basalelement->GetInput(BaseEnum);
+        Input* water_input        = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input);
         if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);                     _assert_(old_wh_input);}
-
-        IssmDouble sediment_storing    = SedimentStoring(basalelement);
 
         /*Transfer related Inputs*/
@@ -326 +327 @@
                 active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
         }
+
 
         /* Start  looping on the number of gaussian points: */
@@ -347 +349 @@
                         /*if EPL is present and active input is there not here*/
                         active_element_input->GetInputValue(&active_element);
-                        if(!active_element){    
+                        if(!active_element){
                                 water_input->GetInputValue(&water_load,gauss);
                                 scalar = Jdet*gauss->weight*(water_load);
@@ -356 +358 @@
                         }
                 }
+
+
                 /*Transient and transfer terms*/
                 if(dt!=0.){
                         old_wh_input    ->GetInputValue(&water_head,gauss);
+                        sediment_storing = SedimentStoring(basalelement,gauss,sed_head_input,base_input);//sed_head_input
                         if(isefficientlayer){
                                 /*Dealing with the sediment part of the transfer term*/
@@ -508 +513 @@
 
 /*Intermediaries*/
-IssmDouble HydrologyDCInefficientAnalysis::SedimentStoring(Element* element){/*{{{*/
+IssmDouble HydrologyDCInefficientAnalysis::SedimentStoring(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input){/*{{{*/
+        IssmDouble sediment_storing;
+        IssmDouble specific_porosity;
+        IssmDouble base_elev,prestep_head,water_sheet;
         IssmDouble rho_freshwater           = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
         IssmDouble g                        = element->GetMaterialParameter(ConstantsGEnum);
@@ -515 +523 @@
         IssmDouble sediment_compressibility = element->GetMaterialParameter(HydrologydcSedimentCompressibilityEnum);
         IssmDouble water_compressibility    = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum);
-        return rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));               
-}/*}}}*/
-
-IssmDouble HydrologyDCInefficientAnalysis::EplSpecificStoring(Element* element){/*{{{*/
-        IssmDouble rho_freshwater        = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
-        IssmDouble g                     = element->GetMaterialParameter(ConstantsGEnum);
-        IssmDouble epl_porosity          = element->GetMaterialParameter(HydrologydcEplPorosityEnum);
-        IssmDouble epl_compressibility   = element->GetMaterialParameter(HydrologydcEplCompressibilityEnum);
-        IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum);
-        return rho_freshwater*g*epl_porosity*(water_compressibility+(epl_compressibility/epl_porosity));                 
+        base_input->GetInputValue(&base_elev,gauss);
+        sed_head_input->GetInputValue(&prestep_head,gauss);
+        water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
+        specific_porosity=sediment_porosity-rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
+        if (water_sheet<=0.9*sediment_thickness){//porosity for unconfined region
+                sediment_storing=sediment_porosity;
+        }
+        else if((water_sheet<sediment_thickness) && (water_sheet>0.9*sediment_thickness)){//continuity ramp
+                sediment_storing=(sediment_thickness-water_sheet)*specific_porosity/(0.1*sediment_thickness)+
+                        rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
+        }
+        else{//storing coefficient for confined
+                sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
+        }
+        return sediment_storing;
+}/*}}}*/
+IssmDouble HydrologyDCInefficientAnalysis::SedimentTransmitivity(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input,Input* SedTrans_input){/*{{{*/
+        IssmDouble sediment_transmitivity;
+        IssmDouble FullLayer_transmitivity;
+        IssmDouble base_elev,prestep_head,water_sheet;
+        IssmDouble sediment_thickness       = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum);
+        base_input->GetInputValue(&base_elev,gauss);
+        sed_head_input->GetInputValue(&prestep_head,gauss);
+        SedTrans_input->GetInputValue(&FullLayer_transmitivity,gauss);
+        water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
+        if (water_sheet<=sediment_thickness){
+                sediment_transmitivity=FullLayer_transmitivity*water_sheet/sediment_thickness;
+        }
+        else{
+                sediment_transmitivity=FullLayer_transmitivity;
+        }
+        return sediment_transmitivity;
 }/*}}}*/
 

issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h

@@ -34 +34 @@
                 /*Intermediaries*/
                 void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-                IssmDouble SedimentStoring(Element* element);
-                IssmDouble EplSpecificStoring(Element* element);
+                IssmDouble SedimentStoring(Element* element, Gauss* gauss, Input* sed_head_input, Input* base_input);
+                IssmDouble SedimentTransmitivity(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input,Input* SedTrans_input);
                 IssmDouble GetHydrologyDCInefficientHmax(Element* element, Gauss* gauss, Input* thickness_input, Input* base_input);
                 void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode);