Changeset 27297

Timestamp:

10/06/22 23:34:44 (2 years ago)

Author:

rueckamp

Message:

NEW: SMB debris model based on Mayer & Licuilli 2021

Location:

issm/trunk-jpl/src

Files:

• c/analyses/SmbAnalysis.cpp (modified) (3 diffs)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (1 diff)
• c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (2 diffs)
• c/shared/Enum/EnumToStringx.cpp (modified) (2 diffs)
• c/shared/Enum/StringToEnumx.cpp (modified) (14 diffs)
• c/shared/io/Marshalling/IoCodeConversions.cpp (modified) (1 diff)
• m/classes/SMBdebrisML.m (added)
• m/classes/initialization.m (modified) (6 diffs)

issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp ¶

@@ -419 +419 @@
                         /*Nothing to add to parameters*/
                         break;
+                case SMBdebrisMLEnum:
+                        /*Nothing to add to parameters*/
+                        break;
                 default:
                         _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
@@ -484 +487 @@
                         break;
                 case SMBarmaEnum:
-         if(VerboseSolution())_printf0_("    call smb arma module\n");
-         Smbarmax(femmodel);
-         break;
+                        if(VerboseSolution())_printf0_("    call smb arma module\n");
+                        Smbarmax(femmodel);
+                        break;
                 case SMBhenningEnum:
                         if(VerboseSolution())_printf0_("  call smb Henning module\n");
@@ -514 +517 @@
                         #endif //_HAVE_SEMIC_
                         break;
+                case SMBdebrisMLEnum:
+                        if(VerboseSolution())_printf0_("        call smb debris Mayer & Liculli module\n");
+                        SmbDebrisMLx(femmodel);
+                        break;
                 default:
                         _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp ¶

@@ -485 +485 @@
                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.rdl",SmbRdlEnum));
                         break;
+                case SMBdebrisMLEnum:
+                        break;
                 default:
                         _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp ¶

@@ -506 +506 @@
 
 }/*}}}*/
+void SmbDebrisMLx(FemModel* femmodel){/*{{{*/
+
+        //      The function is based on:
+        //      Evatt GW, Abrahams ID, Heil M, Mayer C, Kingslake J, Mitchell SL, et al. Glacial melt under a porous debris layer. Journal of Glaciology 61 (2015) 825–836, doi:10.3189/2
+        //      Constants/Values are taken from Mayer, Licciulli (2021): https://www.frontiersin.org/articles/10.3389/feart.2021.710276/full#B7
+        //      function taken from https://github.com/carlolic/DebrisExp/blob/main/USFs/USF_DebrisCoverage.f90
+
+        /*Intermediaries*/
+        // altitude gradients of the crucial parameters (radiation from Marty et al., TaAClimat; 2002)
+        IssmDouble LW=2.9;          // W/m^2 /100m                       2.9
+        IssmDouble SW=1.3;          // W/m^2 /100m                       1.3
+        IssmDouble HumidityG=0;     // % /100m         rough estimate
+        IssmDouble AirTemp=0.7;     // C /100m
+        IssmDouble WindSpeed=0.02;  // m/s /100m       rough estimate    0.2
+
+        // accumulation follows a linear increase above the ELA up to a plateau
+        IssmDouble AccG=0.1;                    // m w.e. /100m
+        IssmDouble AccMax=1.;                    // m w.e.
+        IssmDouble ReferenceElevation=2200.;     // m M&L
+        IssmDouble AblationDays=120.;            //
+
+        IssmDouble In=100.;                 // Wm^-2        incoming long wave
+        IssmDouble Q=500.;                  // Wm^-2        incoming short wave
+        IssmDouble K=0.585;                // Wm^-1K^-1    thermal conductivity          0.585
+        IssmDouble Qm=0.0012;              // kg m^-3      measured humiditiy level
+        IssmDouble Qh=0.006 ;              // kg m^-3      saturated humidity level
+        IssmDouble Tm=2.;                   // C            air temperature
+        IssmDouble Rhoaa=1.22;             // kgm^-3       air densitiy
+        IssmDouble Um=1.5;                 // ms^-1        measured wind speed
+        IssmDouble Xm=1.5;                 // ms^-1        measurement height
+        IssmDouble Xr=0.01;                // ms^-1        surface roughness             0.01
+        IssmDouble Alphad=0.07;            //              debris albedo                 0.07
+        IssmDouble Alphai=0.4;             //              ice ablbedo
+        IssmDouble Ustar=0.16;             // ms^-1        friction velocity             0.16
+        IssmDouble Ca=1000.;                // jkg^-1K^-1   specific heat capacity of air
+        IssmDouble Lm;//=3.34E+05;            // jkg^-1K^-1   latent heat of ice melt
+        IssmDouble Lv=2.50E+06;            // jkg^-1K^-1   latent heat of evaporation
+        IssmDouble Tf=273.;                 // K            water freeezing temperature
+        IssmDouble Eps=0.95;               //              thermal emissivity
+        IssmDouble Rhoi=900.;               // kgm^-3       ice density
+        IssmDouble Sigma=5.67E-08;         // Wm^-2K^-4    Stefan Boltzmann constant
+        IssmDouble Kstar=0.4;              //              von kármán constant
+        IssmDouble Gamma=180.;              // m^-1         wind speed attenuation        234
+        IssmDouble PhiD;//=0.005;              //              debris packing fraction       0.01
+        IssmDouble Humidity=0.2;           //              relative humidity
+
+        IssmDouble smb,yts,z,debris;
+        IssmDouble MassBalanceCmDayDebris,MassBalanceMYearDebris;
+
+        /*Get material parameters and constants */
+        //femmodel->parameters->FindParam(&Rhoi,MaterialsRhoIceEnum); // Note Carlo's model used as  benchmark was run with different densities for debris and FS
+        femmodel->parameters->FindParam(&Lm,MaterialsLatentheatEnum);
+        femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); 
+        PhiD=0.;
+        //if(isdebris) femmodel->parameters->FindParam(&PhiD,DebrisPackingFractionEnum);
+
+        /* Loop over all the elements of this partition */
+        for(Object* & object : femmodel->elements->objects){
+                Element* element=xDynamicCast<Element*>(object);
+
+                /* Allocate all arrays */
+                int         numvertices=element->GetNumberOfVertices();
+                IssmDouble* surfacelist=xNew<IssmDouble>(numvertices);
+                IssmDouble* smb=xNew<IssmDouble>(numvertices);
+                IssmDouble* debriscover=xNew<IssmDouble>(numvertices);
+                element->GetInputListOnVertices(surfacelist,SurfaceEnum);
+
+                /* Get inputs */
+                element->GetInputListOnVertices(debriscover,DebrisThicknessEnum);
+
+                /*Loop over all vertices of element and calculate SMB as function of Debris Cover and z */
+                for(int v=0;v<numvertices;v++){
+
+                        /*Get vertex elevation */
+                        z=surfacelist[v];
+
+                        /* Get debris cover */
+                        debris=debriscover[v];
+
+                        /*IssmDouble dk=1e-5; // TODO make Alphad and Alphai a user input
+                        IssmDouble n=debris/dk;
+                        IssmDouble nmax=1000;
+                        IssmDouble Alphaeff;
+                        if(n>nmax){
+                                Alphaeff=Alphad;
+                        } else {
+                                Alphaeff=Alphai+n*(Alphad-Alphai)/nmax;
+                        }*/
+
+                        // M&L
+                       IssmDouble Alphaeff=Alphad;
+
+                        /* compute smb */
+                        for (int ismb=0;ismb<2;ismb++){
+                                if(ismb==0){
+                                        // calc a reference smb to identify accum and melt region; debris only develops in ablation area
+                                        debris=0.;
+                                }else{
+                                        // only in the meltregime debris develops
+                                        if(-MassBalanceCmDayDebris<0) debris=debriscover[v]; 
+                                }
+                                MassBalanceCmDayDebris=(((In-(z-ReferenceElevation)*LW/100.)-(Eps*Sigma*(Tf*Tf*Tf*Tf))+ 
+                                    (Q+(z-ReferenceElevation)*SW/100.)*(1.-Alphaeff)+ 
+                                    (Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* 
+                                    WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))*(Tm-(z- 
+                                    ReferenceElevation)*AirTemp/100.))/((1-PhiD)*Rhoi*Lm)/(1.+ 
+                                    ((Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* 
+                                    WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/ 
+                                    K*debris)-(Lv*Ustar*Ustar*(Qh-(Qh*(Humidity-(z- 
+                                    ReferenceElevation)*HumidityG/100.)))*(exp(-Gamma*Xr)))/((1.-PhiD)* 
+                                    Rhoi*Lm*Ustar)/((((Um-(z-ReferenceElevation)*WindSpeed/100.) 
+                                    -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)))*100.*24.*60.*60.;
+                        }
+
+                        /* account form ablation days, and convert to m/s */
+                        MassBalanceMYearDebris=-MassBalanceCmDayDebris/100.*AblationDays/yts;
+
+                        /*Update array accordingly*/
+                        smb[v]=MassBalanceMYearDebris;
+                }
+
+                /*Add input to element and Free memory*/
+                element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
+                xDelete<IssmDouble>(surfacelist);
+                xDelete<IssmDouble>(smb);
+                xDelete<IssmDouble>(debriscover);
+        }
+}/*}}}*/
 void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/
 

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h ¶

@@ -23 +23 @@
 void SmbMeltComponentsx(FemModel* femmodel);
 void SmbGradientsComponentsx(FemModel* femmodel);
+void SmbDebrisMLx(FemModel* femmodel);
 /* SEMIC: */
 void SmbSemicx(FemModel* femmodel);

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h ¶

@@ -157 +157 @@
         DamageStressUBoundEnum,
         DebugProfilingEnum,
+        DebrisThicknessEnum,
         DomainDimensionEnum,
         DomainTypeEnum,
@@ -1541 +1542 @@
         SMBarmaEnum,
         SMBcomponentsEnum,
+        SMBdebrisMLEnum,
         SMBd18opddEnum,
         SMBforcingEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp ¶

@@ -165 +165 @@
                 case DamageStressUBoundEnum : return "DamageStressUBound";
                 case DebugProfilingEnum : return "DebugProfiling";
+                case DebrisThicknessEnum : return "DebrisThickness";
                 case DomainDimensionEnum : return "DomainDimension";
                 case DomainTypeEnum : return "DomainType";
@@ -1544 +1545 @@
                 case SMBarmaEnum : return "SMBarma";
                 case SMBcomponentsEnum : return "SMBcomponents";
+                case SMBdebrisMLEnum : return "SMBdebrisML";
                 case SMBd18opddEnum : return "SMBd18opdd";
                 case SMBforcingEnum : return "SMBforcing";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp ¶

@@ -168 +168 @@
               else if (strcmp(name,"DamageStressUBound")==0) return DamageStressUBoundEnum;
               else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum;
+              else if (strcmp(name,"DebrisThickness")==0) return DebrisThicknessEnum;
               else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum;
               else if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
@@ -259 +260 @@
               else if (strcmp(name,"HydrologydcSedimentThickness")==0) return HydrologydcSedimentThicknessEnum;
               else if (strcmp(name,"HydrologyStepAdapt")==0) return HydrologyStepAdaptEnum;
-              else if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"HydrologydcUnconfinedFlag")==0) return HydrologydcUnconfinedFlagEnum;
+              if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum;
+              else if (strcmp(name,"HydrologydcUnconfinedFlag")==0) return HydrologydcUnconfinedFlagEnum;
               else if (strcmp(name,"HydrologyshreveStabilization")==0) return HydrologyshreveStabilizationEnum;
               else if (strcmp(name,"IcecapToEarthComm")==0) return IcecapToEarthCommEnum;
@@ -382 +383 @@
               else if (strcmp(name,"QmuResponsedescriptors")==0) return QmuResponsedescriptorsEnum;
               else if (strcmp(name,"QmuVariableDescriptors")==0) return QmuVariableDescriptorsEnum;
-              else if (strcmp(name,"QmuVariablePartitions")==0) return QmuVariablePartitionsEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"QmuVariablePartitionsNpart")==0) return QmuVariablePartitionsNpartEnum;
+              if (strcmp(name,"QmuVariablePartitions")==0) return QmuVariablePartitionsEnum;
+              else if (strcmp(name,"QmuVariablePartitionsNpart")==0) return QmuVariablePartitionsNpartEnum;
               else if (strcmp(name,"QmuVariablePartitionsNt")==0) return QmuVariablePartitionsNtEnum;
               else if (strcmp(name,"QmuResponsePartitions")==0) return QmuResponsePartitionsEnum;
@@ -505 +506 @@
               else if (strcmp(name,"SmbElevationBins")==0) return SmbElevationBinsEnum;
               else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum;
-              else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"SmbDelta18oSurface")==0) return SmbDelta18oSurfaceEnum;
+              if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum;
+              else if (strcmp(name,"SmbDelta18oSurface")==0) return SmbDelta18oSurfaceEnum;
               else if (strcmp(name,"SmbDenIdx")==0) return SmbDenIdxEnum;
               else if (strcmp(name,"SmbDt")==0) return SmbDtEnum;
@@ -628 +629 @@
               else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
               else if (strcmp(name,"Velocity")==0) return VelocityEnum;
-              else if (strcmp(name,"Xxe")==0) return XxeEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"Yye")==0) return YyeEnum;
+              if (strcmp(name,"Xxe")==0) return XxeEnum;
+              else if (strcmp(name,"Yye")==0) return YyeEnum;
               else if (strcmp(name,"Zze")==0) return ZzeEnum;
               else if (strcmp(name,"Areae")==0) return AreaeEnum;
@@ -751 +752 @@
               else if (strcmp(name,"Domain3D")==0) return Domain3DEnum;
               else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum;
-              else if (strcmp(name,"DrivingStressX")==0) return DrivingStressXEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"DrivingStressY")==0) return DrivingStressYEnum;
+              if (strcmp(name,"DrivingStressX")==0) return DrivingStressXEnum;
+              else if (strcmp(name,"DrivingStressY")==0) return DrivingStressYEnum;
               else if (strcmp(name,"Dummy")==0) return DummyEnum;
               else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum;
@@ -874 +875 @@
               else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum;
               else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
-              else if (strcmp(name,"Misfit")==0) return MisfitEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"MovingFrontalVx")==0) return MovingFrontalVxEnum;
+              if (strcmp(name,"Misfit")==0) return MisfitEnum;
+              else if (strcmp(name,"MovingFrontalVx")==0) return MovingFrontalVxEnum;
               else if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum;
               else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
@@ -997 +998 @@
               else if (strcmp(name,"SmbC")==0) return SmbCEnum;
               else if (strcmp(name,"SmbCcsnowValue")==0) return SmbCcsnowValueEnum;
-              else if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum;
+              if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum;
+              else if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum;
               else if (strcmp(name,"SmbD")==0) return SmbDEnum;
               else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum;
@@ -1120 +1121 @@
               else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
               else if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum;
-              else if (strcmp(name,"Surface")==0) return SurfaceEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
+              if (strcmp(name,"Surface")==0) return SurfaceEnum;
+              else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
               else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
               else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
@@ -1243 +1244 @@
               else if (strcmp(name,"Outputdefinition66")==0) return Outputdefinition66Enum;
               else if (strcmp(name,"Outputdefinition67")==0) return Outputdefinition67Enum;
-              else if (strcmp(name,"Outputdefinition68")==0) return Outputdefinition68Enum;
          else stage=11;
    }
    if(stage==11){
-              if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum;
+              if (strcmp(name,"Outputdefinition68")==0) return Outputdefinition68Enum;
+              else if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum;
               else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum;
               else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum;
@@ -1366 +1367 @@
               else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
               else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
-              else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
          else stage=12;
    }
    if(stage==12){
-              if (strcmp(name,"Element")==0) return ElementEnum;
+              if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
+              else if (strcmp(name,"Element")==0) return ElementEnum;
               else if (strcmp(name,"ElementHook")==0) return ElementHookEnum;
               else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
@@ -1489 +1490 @@
               else if (strcmp(name,"MantlePlumeGeothermalFlux")==0) return MantlePlumeGeothermalFluxEnum;
               else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
-              else if (strcmp(name,"Masscon")==0) return MassconEnum;
          else stage=13;
    }
    if(stage==13){
-              if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
+              if (strcmp(name,"Masscon")==0) return MassconEnum;
+              else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
               else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
               else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
@@ -1580 +1581 @@
               else if (strcmp(name,"SMBarma")==0) return SMBarmaEnum;
               else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
+              else if (strcmp(name,"SMBdebrisML")==0) return SMBdebrisMLEnum;
               else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
               else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
@@ -1611 +1613 @@
               else if (strcmp(name,"Separate")==0) return SeparateEnum;
               else if (strcmp(name,"Seq")==0) return SeqEnum;
-              else if (strcmp(name,"SmbAnalysis")==0) return SmbAnalysisEnum;
-              else if (strcmp(name,"SmbSolution")==0) return SmbSolutionEnum;
          else stage=14;
    }
    if(stage==14){
-              if (strcmp(name,"SmoothAnalysis")==0) return SmoothAnalysisEnum;
+              if (strcmp(name,"SmbAnalysis")==0) return SmbAnalysisEnum;
+              else if (strcmp(name,"SmbSolution")==0) return SmbSolutionEnum;
+              else if (strcmp(name,"SmoothAnalysis")==0) return SmoothAnalysisEnum;
               else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum;
               else if (strcmp(name,"SpatialLinearFloatingMeltRate")==0) return SpatialLinearFloatingMeltRateEnum;

issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp ¶

@@ -240 +240 @@
                 case 10: return SMBpddSicopolisEnum;
                 case 11: return SMBgradientscomponentsEnum;
-                case 12: return SMBsemicEnum;   
+                case 12: return SMBsemicEnum;    
                 case 13: return SMBarmaEnum;
+                case 14: return SMBdebrisMLEnum;
                 default: _error_("Marshalled SMB code \""<<enum_in<<"\" not supported yet");
         }

issm/trunk-jpl/src/m/classes/initialization.m ¶

@@ -25 +25 @@
                 str                 = NaN;
                 sample              = NaN;
+                debris              = NaN;
         end
         methods
@@ -121 +122 @@
                                 if ~isnan(md.initialization.sample)
                                         md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
+                                end
+                        end
+                        if ismember('DebrisAnalysis',analyses),
+                                if ~isnan(md.initialization.debris)
+                                        md = checkfield(md,'fieldname','initialization.debris','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
                                 end
                         end
@@ -145 +151 @@
                         fielddisplay(self,'dsl','Dynamic sea level.');
                         fielddisplay(self,'str','Steric sea level.');
+                        fielddisplay(self,'debris','Surface debris layer [m]');
                 end % }}}
                 function marshall(self,prefix,md,fid) % {{{
@@ -167 +174 @@
                         WriteData(fid,prefix,'object',self,'fieldname','hydraulic_potential','format','DoubleMat','mattype',1);
                         WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1);
+                        WriteData(fid,prefix,'object',self,'fieldname','debris','format','DoubleMat','mattype',1);
 
                         if md.thermal.isenthalpy,
@@ -197 +205 @@
                         self.dsl=project3d(md,'vector',self.dsl,'type','node','layer',1);
                         self.str=project3d(md,'vector',self.str,'type','node','layer',1);
+                        self.str=project3d(md,'vector',self.debris,'type','node','layer',1);
 
                         %Lithostatic pressure by default
@@ -218 +227 @@
                         writejs1Darray(fid,[modelname '.initialization.channel'],self.channelarea);
                         writejs1Darray(fid,[modelname '.initialization.sample'],self.sample);
+                        writejs1Darray(fid,[modelname '.initialization.debris'],self.debris);
 
                 end % }}}