Girinskii potential conversion functions

Language

c++

Purpose

Convert between single layer freshwater/saltwater head and potential. From pages 49, 51, and 100 in Strack, 1989

<math>\Phi=\tfrac{1}{2}k\frac{\rho}{\rho_s-\rho}\left[h-\frac{\rho_s}{\rho}H_{sea}+\frac{\rho_s-\rho}{\rho}H\right]^2</math>

<math>h=\frac{\rho_s}{\rho}H_{sea}-\frac{\rho_s-\rho}{\rho}H+\sqrt{\frac{2}{k}\frac{\rho_s-\rho}{\rho}\Phi}</math>

Where <math>\alpha=\frac{\rho_s}{\rho}</math>, and <math>\rho_s</math> and <math>\rho</math> are the density of saltwater and freshwater, respectively. <math>H_{sea}</math> is the sea level elevation and <math>H</math> is the confined layer thickness.

Author

James R. Craig

#include <math.h>
/********************************************************************************
  IsConfined
   Returns true if local layer condition is confined, false otehrwise
-------------------------------------------------------------------------------*/
inline bool   IsConfined        (const double &potential, 
                                 const double &K, 
                                 const double &H){
  return (potential>=0.5*K*H*H);
}

#include <math.h>
/********************************************************************************
  ConvertToHead
   Converts freshwater potential to LOCAL head (head-base)
   alpha is the ratio of saltwater density to freshwater ~1.035 (1.0 for non-intrusion problems)
   Hsea is the sea level in local head (<=0.0 for non-intrusion problems)
-------------------------------------------------------------------------------*/
inline double ConvertToHead     (const double &potential,
                                 const double &K, 
                                 const double &H,
                                 const double &Hsea,
                                 const double &alpha)
{
  if (K==0){return 0.0;}
  double Hs=max(Hsea,0.0);
  if        (potential>=0.5*K*H*H){                 //confined--------------------------
    if      (potential>=K*H*(alpha*Hs-0.5*H)){      //freshwater head        
      return (potential/(K*H)+ 0.5*H);
    }
    else if (potential<=alpha*K*H*(Hs-0.5*H)){      //saltwater head
      return alpha*Hs-(alpha-1.0)*H;
    }
    else {                                          //mixed head
      return sqrt(2.0*alpha*(potential-alpha*K*H*(Hs-0.5*H))/K*(1-alpha));
    }
  }
  else {                                            //unconfined-------------------------
    if (potential>0.5*K*alpha*alpha*Hs*Hs){         //freshwater head
      return sqrt(2.0*max(potential,0.0)/K);
    }
    else if (potential<0.5*K*alpha*Hs*Hs){          //saltwater head
      return Hs;
    }
    else{                                           //mixed head
      return sqrt(2.0*(alpha-1.0)*(max(potential,0.0)-0.5*K*alpha*Hs*Hs)/(K*alpha))+Hs;
    }
  }
}

#include <math.h>
/********************************************************************************
  ConvertToPotential
   Converts from LOCAL freshwater head to potential
   alpha is the ratio of saltwater density to freshwater ~1.035 (1.0 for non-intrusion problems)
   Hsea is the sea level in local head (<=0.0 for non-intrusion problems)
-------------------------------------------------------------------------------*/
inline double ConvertToPotential(const double &Head, 
                                 const double &K, 
                                 const double &H,
                                 const double &Hsea,
                                 const double &alpha)
{
 double Hs=max(Hsea,0.0);
 if (Head>H){                                    //confined--------------------------
   if (Head>=alpha*Hs){                          //freshwater Head
     return K*H*(Head-0.5*H);
   }
   else{                                         //saltwater Head
     return 0.5*K*(1-alpha)/alpha*pow(Head+(alpha-1)*H-alpha*Hs,2)+alpha*K*H*(Hs-0.5*H); 
   }
 }
 else{                                           //unconfined-------------------------
   if (Head>=alpha*Hs){                          //freshwater Head
     return 0.5*K*Head*Head;
   }
   else{                                         //saltwater Head
     return 0.5*K*(alpha/(alpha-1.0))*pow((Head-Hs),2)+0.5*K*alpha*Hs*Hs;
   }
 }
}