EOS5 / custom capillary-pressure modification: negative pressure and convergence issues when using nonzero entry pressure

Hello everyone,

I am working on a TOUGH3 EOS5 case for a Nuclear waste-cell / FORGE-type model with water + hydrogen and a rectilinear mesh.

I modified the capillary-pressure handling to use a custom negative capillary-pressure convention, and I am trying to keep a nonzero entry pressure (PE) in some materials. I would like help understanding whether the issue is in my constitutive formulation, my EOS modifications, or the way I am handling phase transitions.

What I am doing

I am using EOS5 (water + H2) and a custom PCAP function based on an SGM-type formulation.

In my model:

• PCP(1) is negative

• I am effectively using
  Pl = Pg + Pc

Some materials have a nonzero entry pressure, but Cledz clearly triggers the convergence problem.

• Claym: PE = 2.0E+6

• Bento: PE = 1.33E+6

• Cledz: PE = 0.67E+6

• Dredz: PE = 0.67E+6

The important point is that the issue appears specifically when Cledz uses PE = 0.67E+6.

If I keep Cledz = 0.67E+6, the model develops convergence problems.
If I set Cledz = 0, the model runs well.

So the problem seems to be specifically associated with activating PE = 0.67E+6 in Cledz.

Main symptom

At runtime, I get repeated timestep reductions and then errors like:

NO CONVERGENCE IN SUBROUTINE PP
+++++++++   CANNOT FIND PARAMETERS AT ELEMENT * 554*      XX(M) = -0.60568E+06  0.15715E-06  0.20000E+02

TOUGH ERROR: Error in the Time-dependent bc values.

So the liquid pressure becomes strongly negative because Pc is negative and large in magnitude.

My custom EOS behaviour

In the EOS, the primary variables are updated in the usual way:
 

XX(M)=X(NLOC+M)+XINCR

And I found that XX(1) can become negative during Newton iterations.
 

I tried temporarily protecting pressure with:

XX(1)=MAX(XX(1),1.0D0)

and this keeps the pressure stable enough to continue farther, but I know this is a numerical safeguard and may also affect the solution.

In the two-phase → liquid transition block, I also have logic like:
 

CALL PCAP(SATU,TX,PCP,NPH,NMAT,USRX,K)
pee = -PCP(1)
xx(1)=px-pee
DX(NLOC+1)=XX(1)-X(NLOC+1)

Since PCP(1) is negative, this is equivalent to resetting pressure with px + PCP(1).

My custom PCAP routine

My custom capillary-pressure routine uses a negative capillary-pressure convention. The relevant part is:
 

PC_0=-1.0D0/CP(3,NMAT)
PCE=-CP(6,NMAT)
...
PC=SE**(-1.0D0/AM)-1.0D0
PC=PC**(1.0D0-AM)
PC=PC_0*PC
IF(CP(4,NMAT).NE.0.0D0) PC=MAX(PC,-ABS(CP(4,NMAT)))
PCP(1)=PC

So CP(4) is acting as a soft cap on how negative Pc can become.

Parameter tests I tried

I tested several things:

• setting PE = 0 → model runs well

• keeping PE > 0 → convergence problems return

• lowering CP(4) from 1.0E9 to 1.0E5 helped a bit

• lowering CP(4) too much (for example 1.0E4) caused linear-solver issues

However, even after these changes, I still see PP failures and negative pressures during Newton iterations.

My question

I would appreciate advice on the following:

1. Is this negative-Pc convention fundamentally problematic in EOS5, or can it be made robust?

2. Is the main issue likely to be:

   • my PCAP formulation,

   • the EOS phase-transition reset logic,

   • Or the Newton update of pressure?

3. For people who have implemented custom capillary-pressure laws in TOUGH3:

   • Is it better to limit DX for pressure?

   • protect XX(1),

   • or reformulate the phase-transition pressure reset?

4. Has anyone used a similar nonzero PE formulation successfully with EOS5 / H2–water systems?

Summary

• EOS: EOS5

• Components: water + hydrogen

• The problem appears only when a nonzero entry pressure is used

• PE = 0 runs fine

• Custom PCAP returns negative PCP(1)

• Repeated failures occur in PP, with XX(1) becoming negative

• Sometimes the run later stops with Time-dependent bc values or VISCOH2 range errors

Any suggestions on the best way to keep the physical nonzero PE while avoiding these EOS/convergence failures would be very helpful.

Thank you.