Inverse modeling of HYDRUS-1D.

A discussion forum for Hydrus-1D users.
Post Reply
Mitra
Posts: 16
Joined: Thu Dec 06, 2018 3:44 pm
Location: Switzerland

Inverse modeling of HYDRUS-1D.

Post by Mitra » Tue Jan 14, 2020 3:04 pm

Dear Jirka

I am using HYDRUS-1D to inversely calibrate my hydraulic parameters and the longitudinal diffusivity against soil and tracer data. I have found reasonable numbers with low errors, still, I was interested to check the procedure that I use with you to make sure that I am doing it correctly as in the manual there isn't a step by step description of it. Here is what I do:
1. Entering the observation data in data for the inverse solution (with all manually determined weights set to 1 in this table)
2. starting from a sensible range and initial value for desired parameters to be calibrated.
3. running the inverse modeling, changing the initial value for inverse modeling to a value that falls within the 95% confidence interval, and rerunning the inverse modeling until I find the calibrated parameter and the initial value to stay within the 95% confidence interval.

4. Redoing steps 1-3 with trying different weights (in the weighting of inversion data: standard deviation, mean and no internal weighting) and initial values, or ranges. The results didn't vary much.

To do a check after I found the primary optimal parameters, I tried calibrating parameters only one at a time, and they didn't vary much. Given good statistics from HYDRUS-1D for the inverse modeling, the low variation of calibrated parameters under the above-mentioned conditions, and the good match between observation and my model, I concluded that the calibrated parameters are close to the real values.

I did the aforementioned based on what I understood from the manual, but to be sure that I am on the right track, I am thinking if you would recommend something more? or you would have any advice for me on this matter.

Thanks a lot for your time and care.

Sincerely,
Mitra

Jirka
Posts: 4980
Joined: Sat Mar 16, 2002 3:47 pm
Location: USA
Location: Riverside, CA

Re: Inverse modeling of HYDRUS-1D.

Post by Jirka » Wed Jan 15, 2020 11:30 pm

It seems that you do all the right things.

Obviously, since HYDRUS-1D uses a local gradient optimization method (and not a global optimization method), you could do some more comprehensive analysis. See, for example, some recent papers by Brunetti:

Brunetti, G., J. Šimůnek, and P. Piro, A comprehensive analysis of the variably saturated hydraulic behavior of a green roof in a Mediterranean climate, Vadose Zone Journal, 15(9), pp. 15, doi: 10.2136/vzj2016.04.0032, 2016.

Brunetti, G., J. Šimůnek, and P. Piro, A comprehensive numerical analysis of the hydraulic behavior of a permeable pavement, Journal of Hydrology, 540, 1146-1161, doi: 10.1016/j.jhydrol.2016.07.030, 2016.

Brunetti, G., J. Šimůnek, H. Bogena, R. Baatz, J. A. Huisman, H. Dahlke, and H. Vereecken, On the information content of cosmic-ray neutrons in the inverse estimation of soil hydraulic properties, Vadose Zone Journal, 18, 180123, 24 p., doi: 10.2136/vzj2018.06.0123, 2019.

j.

Mitra
Posts: 16
Joined: Thu Dec 06, 2018 3:44 pm
Location: Switzerland

Re: Inverse modeling of HYDRUS-1D.

Post by Mitra » Tue Jan 21, 2020 3:59 pm

Thanks for your reply.
I have an additional question about the objective function that I couldn't find its answer in the references provided in HYDRUS-1D interface and manual.
If we assume to only have the data on the first term in the objective function introduced in the manual (eq. 9.1). There are 2 weights, v_j, and w_i. v_j weights the different sets of datasets, and w_i that weights the individual points within a dataset. v_j balances the error between datasets and is equivalent to the reciprocal of the variance multiplied by the number of points in each dataset. w_i on the other hand, I assume, is the weight that can be chosen by the user to be equivalent to standard deviation, no internal weighting and mean ratios while v_j formula remains the same. Did I understand this correctly?

Also, I computed the standard deviations for my data set they are respectively 0.36 and 0.006 but the weights in the Fit.out are ~0.214 and ~13.7 (when running the inverse modelling by standard deviation). I was wondering if you could help me to understand how these values of standard deviation have been converted to these weights. The manual weights are 1 for all points.

Jirka
Posts: 4980
Joined: Sat Mar 16, 2002 3:47 pm
Location: USA
Location: Riverside, CA

Re: Inverse modeling of HYDRUS-1D.

Post by Jirka » Tue Jan 21, 2020 4:40 pm

The easiest way to explain is to give you the code (when the standard deviation is used):
else if(iWeight.eq.2) then
if(nQ.gt.0) QAve=QAve/nQ
if(nH.gt.0) hAve=hAve/nh
if(nTh.gt.0) ThAve=ThAve/nTh
if(nV.gt.0) vAve=vAve/nV
if(nC.gt.0) cAve=cAve/nC
if(nK.gt.0) rKAve=rKAve/nK
QSig2=0.
hSig2=0.
ThSig2=0.
vSig2=0.
cSig2=0.
rKSig2=0.
do 6 n=1,NOBb
if(iType(n).eq.0) QSig2=QSig2+(sngl(dabs(FO(n)))-QAve)**2/nQ
if(iType(n).eq.1.or.iType(n).eq.12)
! hSig2=hSig2+(sngl(dabs(FO(n)))-hAve)**2/nh
if(iType(n).eq.2.or.iType(n).eq.5.or.iType(n).eq.13)
! ThSig2=ThSig2+(sngl(FO(n))-ThAve)**2/nTh
if(iType(n).eq.3) vSig2=vSig2+(sngl(dabs(FO(n)))-vAve)**2/nV
if(iType(n).eq.4.or.iType(n).eq.14.or.iType(n).eq.15)
! cSig2=cSig2+(sngl(dabs(FO(n)))-cAve)**2/nc
if(iType(n).eq.6) rKSig2=rKSig2+(sngl(FO(n))-rKAve)**2/nK
6 continue
do 7 n=1,NOBb
if(iType(n).eq.0.and.QSig2.ne.0.)
! WT(n)=WT(n)/(nQ*QSig2)**0.5
if((iType(n).eq.1.or.iType(n).eq.12).and.hSig2.ne.0.)
! WT(n)=WT(n)/(nh*hSig2)**0.5
if((iType(n).eq.2.or.iType(n).eq.5.or.iType(n).eq.13).
! and.ThSig2.ne.0.)
! WT(n)=WT(n)/(nTh*ThSig2)**0.5
if(iType(n).eq.3.and.vSig2.ne.0.)
! WT(n)=WT(n)/(nV*vSig2)**0.5
if((iType(n).eq.4.or.iType(n).eq.14.or.iType(n).eq.15).
! and.cSig2.ne.0.)
! WT(n)=WT(n)/(nC*cSig2)**0.5
if(iType(n).eq.6.and.rKSig2.ne.0.)
! WT(n)=WT(n)/(nK*rKSig2)**0.5
7 continue
end if

Mitra
Posts: 16
Joined: Thu Dec 06, 2018 3:44 pm
Location: Switzerland

Re: Inverse modeling of HYDRUS-1D.

Post by Mitra » Tue Jan 21, 2020 5:03 pm

Thanks a lot for sharing the code, I will highly appreciate it if you could confirm whether my understanding on v_j is correct?
For example, if I have 2 datasets one with a mean of 0.3 and 10 points and the other with a mean of 0.1 and 20 points, they are weighted by the values of (1/3) and (1/2) respectively, right? and in addition to this weight, they are weighted by the standard deviation or mean ratio or none at all.

Jirka
Posts: 4980
Joined: Sat Mar 16, 2002 3:47 pm
Location: USA
Location: Riverside, CA

Re: Inverse modeling of HYDRUS-1D.

Post by Jirka » Tue Jan 21, 2020 6:55 pm

here is the entire routine. Figure it out. I do not have time for that any more. J.

*************************************************************************

subroutine Weight(HO,FO,WT,iType,iPos,NOBb,NOB,NOB1,NOB2,iWeight,
! Type,iConcType,ierr)

double precision FO(NOBb),WT(NOBb)
dimension HO(NOBb),iType(NOBb),iPos(NOBb)
character TYPE*9

write(6,1001,err=901) Type,Type
n1=0
n2=0
nQ=0
nV=0
nH=0
nTh=0
nC=0
nK=0
QAve=0.0
vAve=0.0
hAve=0.0
cAve=0.0
ThAve=0.0
rKAve=0.0
read(5,*,err=902)
do 4 n=1,NOBb
read(5,*,err=902) HO(n),FOS,iType(n),iPos(n),WTS
if(WTS.eq.0.0) WTS=1.0
FO(n)=dble(FOS)
c logarithm transformation of pressure heads
ihType=0 ! standard approach
c ihType=1 ! log transformation of the pressure heads
if(ihType.eq.1.and.iType(n).eq.1) then
if(abs(FO(n)).gt.1.d-15) then
FO(n)=dlog10(abs(FO(n)))
else
FO(n)=-15.
end if
end if
c logarithm transform of resident and total concentrations
if(iConcType.eq.1.and.iType(n).eq.4.and.iPos(n).ne.0) then
if(FO(n).gt.1.d-15) then
FO(n)=dlog10(FO(n))
else
FO(n)=-15.
end if
end if
if(iConcType.eq.1.and.iType(n).eq.15) then
if(FO(n).gt.1.d-15) then
FO(n)=dlog10(FO(n))
else
FO(n)=-15.
end if
end if
WT(n)=dble(WTS)
if(iType(n).le.4.or. iType(n).ge.12) n1=n1+1
if(iType(n).ge.7.and.iType(n).le.11) n2=n2+1
if(iType(n).eq.0) then
nQ=nQ+1
QAve=QAve+abs(sngl(FO(n)))
else if (iType(n).eq.1.or.iType(n).eq.12) then
nH=nH+1
hAve=hAve+abs(sngl(FO(n)))
else if (iType(n).eq.2.or.iType(n).eq.5.or.iType(n).eq.13) then
nTh=nTh+1
ThAve=ThAve+abs(sngl(FO(n)))
else if (iType(n).eq.3) then
nV=nV+1
vAve=vAve+abs(sngl(FO(n)))
else if (iType(n).eq.4.or.iType(n).eq.14.or.iType(n).eq.15) then
nC=nC+1
cAve=cAve+abs(sngl(FO(n)))
else if (iType(n).eq.6) then
nK=nK+1
rKAve=rKAve+abs(sngl(FO(n)))
end if
4 continue
NOB=n1
NOB1=n2
NOB2=NOBb-NOB-NOB1

* ----- Default weighting for observation data ----
if(iWeight.eq.1) then
Summ=1.
if(nQ+nH+nTh+nK+nV+nC.gt.0)
! Summ=(QAve+hAve+ThAve+rKAve+vAve+cAve)/(nQ+nH+nTh+nK+nV+nC)
if(nQ.gt.0) then
QAve=QAve/nQ
Summ=amin1(Summ,QAve)
end if
if(nH.gt.0) then
hAve=hAve/nh
Summ=amin1(Summ,hAve)
end if
if(nTh.gt.0) then
ThAve=ThAve/nTh
Summ=amin1(Summ,ThAve)
end if
if(nV.gt.0) then
vAve=vAve/nV
Summ=amin1(Summ,vAve)
end if
if(nC.gt.0) then
cAve=cAve/nC
Summ=amin1(Summ,cAve)
end if
if(nK.gt.0) then
rKAve=rKAve/nK
Summ=amin1(Summ,rKAve)
end if
do 5 n=1,NOBb
if(iType(n).eq.0.and.QAve.ne.0.) WT(n)=WT(n)/QAve*Summ
if((iType(n).eq.1.or.iType(n).eq.12).and.hAve.ne.0.)
! WT(n)=WT(n)/hAve*Summ
if((iType(n).eq.2.or.iType(n).eq.5.or.iType(n).eq.13).and.
! ThAve.ne.0.) WT(n)=WT(n)/ThAve*Summ
if(iType(n).eq.3.and.vAve.ne.0.) WT(n)=WT(n)/vAve*Summ
if((iType(n).eq.4.or.iType(n).eq.14.or.iType(n).eq.15).and.
! cAve.ne.0.) WT(n)=WT(n)/cAve*Summ
if(iType(n).eq.6.and.rKAve.ne.0.) WT(n)=WT(n)/rKAve*Summ
5 continue
else if(iWeight.eq.2) then
if(nQ.gt.0) QAve=QAve/nQ
if(nH.gt.0) hAve=hAve/nh
if(nTh.gt.0) ThAve=ThAve/nTh
if(nV.gt.0) vAve=vAve/nV
if(nC.gt.0) cAve=cAve/nC
if(nK.gt.0) rKAve=rKAve/nK
QSig2=0.
hSig2=0.
ThSig2=0.
vSig2=0.
cSig2=0.
rKSig2=0.
do 6 n=1,NOBb
if(iType(n).eq.0) QSig2=QSig2+(sngl(dabs(FO(n)))-QAve)**2/nQ
if(iType(n).eq.1.or.iType(n).eq.12)
! hSig2=hSig2+(sngl(dabs(FO(n)))-hAve)**2/nh
if(iType(n).eq.2.or.iType(n).eq.5.or.iType(n).eq.13)
! ThSig2=ThSig2+(sngl(FO(n))-ThAve)**2/nTh
if(iType(n).eq.3) vSig2=vSig2+(sngl(dabs(FO(n)))-vAve)**2/nV
if(iType(n).eq.4.or.iType(n).eq.14.or.iType(n).eq.15)
! cSig2=cSig2+(sngl(dabs(FO(n)))-cAve)**2/nc
if(iType(n).eq.6) rKSig2=rKSig2+(sngl(FO(n))-rKAve)**2/nK
6 continue
do 7 n=1,NOBb
if(iType(n).eq.0.and.QSig2.ne.0.)
! WT(n)=WT(n)/(nQ*QSig2)**0.5
if((iType(n).eq.1.or.iType(n).eq.12).and.hSig2.ne.0.)
! WT(n)=WT(n)/(nh*hSig2)**0.5
if((iType(n).eq.2.or.iType(n).eq.5.or.iType(n).eq.13).
! and.ThSig2.ne.0.)
! WT(n)=WT(n)/(nTh*ThSig2)**0.5
if(iType(n).eq.3.and.vSig2.ne.0.)
! WT(n)=WT(n)/(nV*vSig2)**0.5
if((iType(n).eq.4.or.iType(n).eq.14.or.iType(n).eq.15).
! and.cSig2.ne.0.)
! WT(n)=WT(n)/(nC*cSig2)**0.5
if(iType(n).eq.6.and.rKSig2.ne.0.)
! WT(n)=WT(n)/(nK*rKSig2)**0.5
7 continue
end if
call SortData(NOBb,HO,FO,WT,iPos,iType)
do 8 n=1,NOBb
if(abs(sngl(FO(n))).lt.0.01.or.abs(sngl(FO(n))).gt.99999.) then
if(abs(HO(n)).lt.0.01.or.abs(HO(n)).gt.99999.) then
write(6,1002,err=901) n,HO(n),sngl(FO(n)),iType(n),iPos(n),
! sngl(WT(n))
else
write(6,1003,err=901) n,HO(n),sngl(FO(n)),iType(n),iPos(n),
! sngl(WT(n))
end if
else
if(abs(HO(n)).lt.0.01.or.abs(HO(n)).gt.99999.) then
write(6,1004,err=901) n,HO(n),sngl(FO(n)),iType(n),iPos(n),
! sngl(WT(n))
else
write(6,1005,err=901) n,HO(n),sngl(FO(n)),iType(n),iPos(n),
! sngl(WT(n))
end if
end if
8 continue
return

901 ierr=1
return
902 ierr=2
return

1001 format(//1X,'Observed ',A/1X,18(1H=)/4X,'Obs',5X,'Time',4X,A,4X,
!'Type',2X,'Position'3X,'Weight')
1002 format(1X,I5,2E11.3,5X,I3,I6,F14.6)
1003 format(1X,I5,F11.3,E11.3,5X,I3,I6,F14.6)
1004 format(1X,I5,e11.3,F11.3,5X,I3,I6,F14.6)
1005 format(1X,I5,2F11.3,5X,I3,I6,F14.6)
end

*************************************************************************

Post Reply