How to implement sorption in unsaturated domain with OGS-5?

Hi everyone,

is anyone familiar with implementing a sorption isotherm under unsaturated conditions? I got some weird results here (amount of sorbed species increaseing more than amount of solved species decreases, causing increase in total amount) and looking into the code in rfmat_cp.cpp it seems like saturation is not considered in computation of the retardation factor.

I’m not familiar with C++ but how I understood the code, the retardation factor seems to be calculated as

retard = 1. + (1. - porosity) * density_rock * isotherm /porosity;

(rfmat_cp.cpp, line 1564 ff) which assumes full saturation by using porosity instead of volumetric water content or did I miss something there? Porosity is said to be a constant value ($POROSITY 1 0.392) and I use a Freundlich isotherm ($ISOTHERM 4 KF nF).
How can I correctly account for sorption under unsaturated conditions?

Cheers, Erik

I think there are no unconfined benchmark examples included but to my understanding you could consider to use a different porosity model.

Now you use $POROSITY 1 0.392, i.e. porosity model = 1. In other words the porosity is constant. Instead you could use e.g. a curve input with case 0. Case 13 look also interesting and it my takes saturation already into account.

Here a list of the possible methods for porosity models:

switch (porosity_model)
	case 0: // n = f(x)
		porosity = GetCurveValue(fct_number, 0, primary_variable[0], &gueltig);
	case 1: // n = const
		porosity = porosity_model_values[0];
	case 2: // n = f(sigma_eff), Stress dependance
		porosity = PorosityEffectiveStress(number, primary_variable[0]);
	case 3: // n = f(S), Free chemical swelling
		    = PorosityVolumetricFreeSwellingConstantIonicstrength(number, primary_variable[0], primary_variable[1]);
	case 4: // n = f(S), Constrained chemical swelling
		porosity = PorosityEffectiveConstrainedSwellingConstantIonicStrength(number, primary_variable[0],
		                                                                     primary_variable[1], &porosity_sw);
	case 5: // n = f(S), Free chemical swelling, I const
		porosity = PorosityVolumetricFreeSwelling(number, primary_variable[0], primary_variable[1]);
	case 6: // n = f(S), Constrained chemical swelling, I const
		    = PorosityEffectiveConstrainedSwelling(number, primary_variable[0], primary_variable[1], &porosity_sw);
	case 7: // n = f(mean stress) WW
		gval = ele_value_dm[number];
		primary_variable[0] = -gval->MeanStress(assem->gp) / 3.0;
		porosity = GetCurveValue(porosity_curve, 0, primary_variable[0], &gueltig);
	case 10:
		/* porosity change through dissolution/precipitation */
		porosity = PorosityVolumetricChemicalReaction(number);
	case 11: // n = temp const, but spatially distributed CB
		// porosity = porosity_model_values[0];
		porosity = _mesh->ele_vector[number]->mat_vector(por_index);
	case 12: // n = n0 + vol_strain, WX: 03.2011
		porosity = PorosityVolStrain(number, porosity_model_values[0], assem);
	case 13:
		CRFProcess* m_pcs_flow = PCSGetFlow();
		const int idx_n = m_pcs_flow->GetElementValueIndex("POROSITY");

		// porosity change through dissolution/precipitation
		// Here, you should access porosity from the element value vector of the flow process
		// so you have to get the index of porosity above, if porosity model = 13
		porosity = m_pcs_flow->GetElementValue(number, idx_n + 1);

Thanks for the answer!

I’m not sure if I’m on the right track, but with assignig a curve in case 0 I would give the water content values in place of the porosity to correctly account for the undersaturation?
Assigning it that way would make porosity (resp. water cont.) time dependent, right? But for the transient flow simulation I want to implement I wouldn’t know the temporal course, plus I would need the assigned values to vary in space…

Case 13 looks interesting, how is it assigned? Simply
(it doesn’t seem to call porosity_model_values[0])

or rather
13 0.392
Do you know how saturation would be considered in this case?

Best, Erik

examples are e.g. in

I never used the method by myself. That’s why I said: it could be possible that it uses saturation but I didn’t look up the method in the code yet.

I think the curve is related to your value of your primary variable :wink:

Thanks for the tips, I’ll give it a shot with the $POROSITY case13