Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 25940)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 25941)
@@ -112,5 +112,5 @@
 	}
 	iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
-	//iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
+	iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
 
 
@@ -153,5 +153,7 @@
 void StressbalanceVerticalAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
 
-	/*No specific parameters*/
+	/* specific parameters*/
+	parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.isvzsolution",StressbalanceIsVzSolutionEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.isvelbc",StressbalanceIsVelBCEnum));
 
 }/*}}}*/
@@ -516,4 +518,7 @@
 	IssmDouble   rho_ice,g;
 
+	int          isvzsolution;
+	IssmDouble   isvelbc, rheology_n, eta;
+
 	/*Get the approximation and do nothing if the element in FS or None*/
 	element->GetInputValue(&approximation,ApproximationEnum);
@@ -521,4 +526,8 @@
 		return;
 	}
+
+	/*Get parameters*/
+	element->FindParam(&isvzsolution,StressbalanceIsVzSolutionEnum);
+	element->FindParam(&isvelbc,StressbalanceIsVelBCEnum);
 
 	/*Get dof list and vertices coordinates: */
@@ -535,4 +544,7 @@
 	IssmDouble*  pressure  = xNew<IssmDouble>(numnodes);
 	IssmDouble*  surface   = xNew<IssmDouble>(numnodes);
+	IssmDouble*  thickness = xNew<IssmDouble>(numnodes);
+	IssmDouble*  base	     = xNew<IssmDouble>(numnodes);
+	IssmDouble*  smb       = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector vz: */
@@ -575,4 +587,29 @@
 		}
 	}
+
+	/*Set analytical vertical velocity field at slow flow region.*/
+	if (isvzsolution == 1){
+		rheology_n=element->material->GetN();
+		element->GetInputListOnNodes(&thickness[0],ThicknessEnum,0.);
+		element->GetInputListOnNodes(&base[0],BaseEnum,0.);
+		element->GetInputListOnNodes(&smb[0],SmbMassBalanceEnum,0.);
+		for(i=0;i<numnodes;i++){
+			eta = (base[i]+thickness[i]-xyz_list[i*3+2])/thickness[i];
+			vz[i] = -(pow(eta,rheology_n+2.0)-1-(rheology_n+2.0)*(eta-1.))/(rheology_n+1.)*smb[i];
+		}	
+	}
+	else if (isvzsolution == 2) {
+		rheology_n=element->material->GetN();
+		element->GetInputListOnNodes(&thickness[0],ThicknessEnum,0.);
+		element->GetInputListOnNodes(&base[0],BaseEnum,0.);
+		element->GetInputListOnNodes(&smb[0],SmbMassBalanceEnum,0.);
+		for(i=0;i<numnodes;i++){
+			if ( vx[i]*vx[i]+vy[i]*vy[i] < isvelbc*isvelbc ){
+				eta = (base[i]+thickness[i]-xyz_list[i*3+2])/thickness[i];
+				vz[i] = -(pow(eta,rheology_n+2.0)-1-(rheology_n+2.0)*(eta-1.))/(rheology_n+1.)*smb[i];
+			}
+		}	
+	}
+
 
 	/*Now Compute vel*/
@@ -600,4 +637,7 @@
 
 	/*Free ressources:*/
+	xDelete<IssmDouble>(smb);
+	xDelete<IssmDouble>(base);
+	xDelete<IssmDouble>(thickness);
 	xDelete<IssmDouble>(surface);
 	xDelete<IssmDouble>(pressure);
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 25940)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 25941)
@@ -441,4 +441,6 @@
 	StressbalanceRiftPenaltyThresholdEnum,
 	StressbalanceShelfDampeningEnum,
+	StressbalanceIsVzSolutionEnum,
+	StressbalanceIsVelBCEnum,
 	ThermalIsdrainicecolumnEnum,
 	ThermalIsdynamicbasalspcEnum,
Index: /issm/trunk-jpl/src/m/classes/stressbalance.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/stressbalance.m	(revision 25940)
+++ /issm/trunk-jpl/src/m/classes/stressbalance.m	(revision 25941)
@@ -23,4 +23,6 @@
 		loadingforce           = NaN;
 		requested_outputs      = {};
+		isvzsolution           = NaN;
+		isvelbc                = NaN;	
 	end
 	methods
@@ -75,4 +77,8 @@
 			 self.rift_penalty_lock=10;
 
+			 %Set vertical velocity option
+			 self.isvzsolution = 0;
+			 self.isvelbc      = 10.;
+
 			 %output default:
 			 self.requested_outputs={'default'};
@@ -96,4 +102,6 @@
 			md = checkfield(md,'fieldname','stressbalance.referential','size',[md.mesh.numberofvertices 6]);
 			md = checkfield(md,'fieldname','stressbalance.loadingforce','size',[md.mesh.numberofvertices 3]);
+			md = checkfield(md,'fieldname','stressbalance.isvzsolution','numel',[1],'values',[0 1 2]);
+			md = checkfield(md,'fieldname','stressbalance.isvelbc','numel',[1],'>=',0.0);
 			md = checkfield(md,'fieldname','stressbalance.requested_outputs','stringrow',1);
 			if ~any(isnan(md.stressbalance.vertex_pairing)),
@@ -158,4 +166,8 @@
 			fielddisplay(self,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset');
 			fielddisplay(self,'vertex_pairing','pairs of vertices that are penalized');
+
+			disp(sprintf('\n      %s','Estimation of vertical velocity field options:'));
+			fielddisplay(self,'isvzsolution','0: incompressible assumption(IA), 1: internal deformation(ID), 2: IA+ID');
+			fielddisplay(self,'isvelbc','If turned on isvzsolution as 2, region (vel < isvelbc) calculates for ID.');
 
 			disp(sprintf('\n      %s','Other:'));
@@ -187,4 +199,6 @@
 			WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','rift_penalty_threshold','format','Integer');
 			WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','referential','format','DoubleMat','mattype',1);
+			WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','isvzsolution','format','Integer','numel',1,'values',[0 1 2]);
+			WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','isvelbc','format','Double','numel',1,'scale',1./yts);
 
 			if size(self.loadingforce,2)==3,