How to implement mixed Ditrichlet and Neumann boundary conditions into pdepe-solver, for solving richards equation with seepage face?

Hi, I have the attached code for solving the richards equation using the pdepe-solver (adapted from link). I am struggling with the boundary conditions. I want to model water content in a soil column above a drain, where (when the bottom soil layer is fully saturated) water can pass downwards through the bottom boundary, but never upwards. The linked document recommends using a conditional bottom boundary (seepage face), with Dirichlet condition when the bottom layer is saturated (hr=0), and Neuman condition (no-pass) when the bottom layer is unsaturated (hr<0). Do you have any idea on how I could implement these conditions into the pdepe-solver? As of now the model uses only Dirichlet conditions for the bottom boundary with a fixed value of hr=0.

L       = 90;                   % soil depth [cm]
dz      = 3;                    % distance between mesh points [cm]
mesh    = L/dz+1;               % number of mesh points
z       = linspace(0,L3,mesh);  % spatial variable [cm]
T       = 72;                   % total time [h]
t       = linspace(0,T,T+1);    % time varable [h]
v_ir    = [1.5,zeros(1,T)];     % irrigation velocity [cm/h]
v_et    = 0.015*ones(1,T+1);    % evapotranspiration velocity [cm/h]
h0      = linspace(-50,0,mesh); % initial pressure head [cm]
theta_r = 0.227;                % residual water content [cm/cm]
theta_s = 0.56;                 % saturated water content [cm/cm]
alpha   = 0.11148;              % van Genuchten parameter [1/cm]
lambda  = 0.0001;               % van Genuchten parameter [-]
n       = 1.55938;              % van Genuchten parameter [-]
ks      = 103.46;               % conductivity at saturation [cm/h]


% PDE solver
options=odeset('RelTol',1e-4,...
               'AbsTol',1e-4,...
               'NormControl','off',...
               'InitialStep',1e-7);
h = pdepe(0,@unsatpde,@unsatic,...  % pressure head [cm]
          @unsatbc,z,t,options,...
          v_ir,v_et,theta_r,theta_s,alpha,lambda,n,ks,L,T,h0);
      

% Partial derivative
function [c,f,s] = unsatpde(z,t,h,DhDz,v_ir,v_et,theta_r,theta_s,alpha,lambda,n,ks,L,T,h0)
[theta,k,c] = sedprop(h,theta_r,theta_s,alpha,lambda,n,ks);
f = k.*DhDz-k;  % flux term [cm/h]
s = 0;          % source/sink term [1/h]
if z <= 30      % evapotranspiration evenly distributed over first 30cm
    s = -interp1(linspace(0,T,T+1),v_et,t)/30;
end
end

% initial condition
function h0z = unsatic(z,v_ir,v_et,theta_r,theta_s,alpha,lambda,n,ks,L,T,h0)
mesh    = length(h0);
h0z     = interp1(linspace(0,L,mesh),h0,z);
end

% boundary condition
function [pl,ql,pr,qr] = unsatbc(zl,hl,zr,hr,t,v_ir,v_et,theta_r,theta_s,alpha,lambda,n,ks,L,T,h0)
pl = interp1(linspace(0,T,T+1),v_ir,t,'previous');            
ql = 1; 
pr = hr;
qr = 0;
end

% soil hydraulic properties
function [theta,k,c] = sedprop(h,theta_r,theta_s,alpha,lambda,n,ks)
m = 1-1/n;
if h >= 0   % saturated
    c = 1e-20;
    k = ks;
    theta = theta_s;
else        % unsaturated
    S = (1+(-alpha*h)^n)^-m;                % saturation [-]
    theta = theta_r+(theta_s-theta_r)*S;    % soil water content [-]
    c = ((theta_s-theta_r)*n*m*alpha*(-alpha*h)^(n-1))/...  % soil retention curve [1/cm]
        ((1+(-alpha*h)^n)^(m+1))+1.e-20;
    k = ks*S^lambda*(1-(1-S^(1/m))^m)^2;    % conductivity [-]
end
end