performing wetting condition on Twophaseflow engine

Hi every body,

I am going to simulate wetting condition after drying with employing Twophaseflow engine. I mention my code in below, at the
stage two when the capillary pressure (PC) decrease slightly it is expected that the degree of saturation increase slowly during the wetting, but degree of saturation was not changed and it is constant at the two stage !

would you please, give me the best way to simulate wetting condition?
what should I do to simulate wetting exactly after the drying according to the below code ?

#######################################################
## Stage one (Drainage) ###
#######################################################
triax.stressMask=7
triax.goal1=triax.goal3=confiningS
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.wall_bottom_activated=0
recorder.dead=0
VTK.dead=1

waterContentMatrix=[];suctionMatrix=[]

print ('voidratio=',(triax.porosity)/(1-triax.porosity))

unsat=TwoPhaseFlowEngine()
unsat.entryPressureMethod=1
meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius) / 2.

print ('meanDiameter=',meanDiameter)

##set boundary conditions, the drainage is controlled by decreasing W-phase pressure and keeping NW-phase pressure constant
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondValue=[0,0,-1e8,0,0,0]
unsat.isPhaseTrapped=True #the W-phase can be disconnected from its reservoir
unsat.initialization()
unsat.surfaceTension =10

file=open('pcSwStrain'+key+'.txt',"w")
for pg in arange(1e-5,25,5):
  unsat.bndCondValue=[0,0,(-1.0)*pg*unsat.surfaceTension/meanDiameter,0,0,0]

  unsat.invasion()

  VolumetricWaterContent=unsat.getSaturation(False)*(triax.porosity)
  waterContentMatrix.append(VolumetricWaterContent)
  suctionMatrix.append(-unsat.bndCondValue[2])

  unsat.computeCapillaryForce()
  for b in O.bodies:
    O.forces.setPermF(b.id, unsat.fluidForce(b.id))
  while 1:
    O.run(1000,True)
    unb=unbalancedForce()
    print ('e=', triax.porosity/(1-triax.porosity), 'suction=' , (-unsat.bndCondValue[2]*1e-3), 'w=', ((unsat.getSaturation(False)*(triax.porosity/(1-triax.porosity)))/2.73), 's22=' ,(-triax.stress(triax.wall_top_id)[1]), 'e11=' ,(-triax.strain[1]), 's33=' , (-triax.stress(triax.wall_front_id)[2]), 'Sr=', (unsat.getSaturation(False)),'v=' ,VolumetricWaterContent)
    if unb<0.1:
      break

print("stage two")
##########################
#stage two (wetting)
##########################
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondValue=[0,0,-1e8,0,0,0]

unsat.initialization()
unsat.surfaceTension = 10
unsat.isPhaseTrapped=True
unsat.isSideBoundaryIncluded=True
unsat.isInvadeBoundary=True
for pg2 in arange(20,1e-5,-5):
  unsat.bndCondValue=[0,0,(-1.0)*pg2*unsat.surfaceTension/meanDiameter,0,0,0]
  unsat.invasion()
  #Suction.append((-unsat.bndCondValue[2])*1e-3)
  unsat.computeCapillaryForce()
  for b in O.bodies:
    O.forces.setPermF(b.id, unsat.fluidForce(b.id))
  while 1:
    O.run(1000,True)
    unb=unbalancedForce()
    if unb<0.25:
      break
  print ('PC:',((-unsat.bndCondValue[2])*1e-3), 'Sr:',(unsat.getSaturation(False)),'e:',(triax.porosity)/(1-triax.porosity), 'e33=', (triax.strain[2]), 'e22=', (-triax.strain[1]), 's22=', (-triax.stress(3)[1]*1e-3))