Skip to content

Waq/bugfix/delwaq 1232 improve settling deposition #561

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
arjenmarkusdeltaresnl wants to merge 6 commits into
base: main
Choose a base branch
from
Open

Waq/bugfix/delwaq 1232 improve settling deposition #561

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
3f4dcd9
DELWAQ-1232: A very first step towards solving the deposition issue
arjenmarkusdeltaresnl Jan 16, 2026
8dd119f
DELWAQ-1232: Adjust the ad hoc code for a column model
arjenmarkusdeltaresnl Jan 19, 2026
9914093
Merge branch 'main' into waq/bugfix/DELWAQ-1232-improve-settling-depo…
arjenmarkusdeltaresnl Jan 28, 2026
4dcd750
Generalise the robust implementation of sedimentaiton and deposition
arjenmarkusdeltaresnl Jan 28, 2026
31be803
DELWAQ-1232: correct calculation of maximum allowable flux
arjenmarkusdeltaresnl Jan 30, 2026
9a65dd9
DELWAQ-1232: Implement a more sophisticated estimate of the maximum m…
arjenmarkusdeltaresnl Feb 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 12 additions & 3 deletions src/engines_gpl/waq/resources/process_lib/csvFiles/inputs.csv
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3505,7 +3505,10 @@ inpupr,inpuit,inpunm,inpude,inpudo,inpusx
"Sed_IM1 ","FrTIMS2 ", 10,"Y","x",1
"Sed_IM1 ","FrTIMS2Max", 11,"Y","x",1
"Sed_IM1 ","PsedminIM1", 12,"Y","x",1
"Sed_IM1 ","VSedIM1 ", 13,"Y"," ",0
"Sed_IM1 ","Volume ", 13,"N","x",1
"Sed_IM1 ","WSSedim ", 14,"N"," ",1
"Sed_IM1 ","VSedIM1 ", 15,"Y"," ",0
"Sed_IM1 ","Flow ", 16,"Y","x",0
"Sed_IM2 ","IM2 ", 1,"N","x",1
"Sed_IM2 ","ZSedIM2 ", 2,"Y","x",1
"Sed_IM2 ","VSedIM2 ", 3,"Y","x",1
Expand All @@ -3518,7 +3521,10 @@ inpupr,inpuit,inpunm,inpude,inpudo,inpusx
"Sed_IM2 ","FrTIMS2 ", 10,"Y","x",1
"Sed_IM2 ","FrTIMS2Max", 11,"Y","x",1
"Sed_IM2 ","PsedminIM2", 12,"Y","x",1
"Sed_IM2 ","VSedIM2 ", 13,"Y"," ",0
"Sed_IM2 ","Volume ", 13,"N","x",1
"Sed_IM2 ","WSSedim ", 14,"N"," ",1
"Sed_IM2 ","VSedIM2 ", 15,"Y"," ",0
"Sed_IM2 ","Flow ", 16,"Y","x",0
"Sed_IM3 ","IM3 ", 1,"N","x",1
"Sed_IM3 ","ZSedIM3 ", 2,"Y","x",1
"Sed_IM3 ","VSedIM3 ", 3,"Y","x",1
Expand All @@ -3531,7 +3537,10 @@ inpupr,inpuit,inpunm,inpude,inpudo,inpusx
"Sed_IM3 ","FrTIMS2 ", 10,"Y","x",1
"Sed_IM3 ","FrTIMS2Max", 11,"Y","x",1
"Sed_IM3 ","PsedminIM3", 12,"Y","x",1
"Sed_IM3 ","VSedIM3 ", 13,"Y"," ",0
"Sed_IM3 ","Volume ", 13,"N","x",1
"Sed_IM3 ","WSSedim ", 14,"N"," ",1
"Sed_IM3 ","VSedIM3 ", 15,"Y"," ",0
"Sed_IM3 ","Flow ", 16,"Y","x",0
"SedDiat ","Diat ", 1,"N","x",1
"SedDiat ","ZSedDiat ", 2,"Y","x",1
"SedDiat ","VSedDiat ", 3,"Y","x",1
Expand Down
1 change: 1 addition & 0 deletions src/engines_gpl/waq/resources/process_lib/csvFiles/items.csv
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14360,3 +14360,4 @@ itemid,itemse,itemex,itemde,itemun,itemnm,itemag,itemda,itemwk,itemgr
"dRS1IM2Buf","x"," ", -999.000 ,"(g/m3/d) ","resuspension flux IM2 from layer S1 ","volume ","volume "," "," "
"dRS1IM3Buf","x"," ", -999.000 ,"(g/m3/d) ","resuspension flux IM3 from layer S1 ","volume ","volume "," "," "
"BloomStep ","x"," ", 1.00000 ,"(d) ","Time step for BLOOM ","volume ","volume "," "," "
"WSSedim ","x"," ", -999.000 ,"(-) ","workspace array for sedimentation ","volume ","volume "," "," "
9 changes: 6 additions & 3 deletions src/engines_gpl/waq/resources/process_lib/csvFiles/outputs.csv
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1061,15 +1061,18 @@ outppr,outpit,outpnm,outpdo,outpsx
"Sed_IM1 ","PSedIM1 ", 1,"x",1
"Sed_IM1 ","fSedIM1 ", 2,"x",1
"Sed_IM1 ","fSedIM1S2 ", 3," ",1
"Sed_IM1 ","VxSedIM1 ", 4,"x",0
"Sed_IM1 ","WSSedim ", 4," ",1
"Sed_IM1 ","VxSedIM1 ", 5,"x",0
"Sed_IM2 ","PSedIM2 ", 1,"x",1
"Sed_IM2 ","fSedIM2 ", 2,"x",1
"Sed_IM2 ","fSedIM2S2 ", 3," ",1
"Sed_IM2 ","VxSedIM2 ", 4,"x",0
"Sed_IM2 ","WSSedim ", 4," ",1
"Sed_IM2 ","VxSedIM2 ", 5,"x",0
"Sed_IM3 ","PSedIM3 ", 1,"x",1
"Sed_IM3 ","fSedIM3 ", 2,"x",1
"Sed_IM3 ","fSedIM3S2 ", 3," ",1
"Sed_IM3 ","VxSedIM3 ", 4,"x",0
"Sed_IM3 ","WSSedim ", 4," ",1
"Sed_IM3 ","VxSedIM3 ", 5,"x",0
"SedDiat ","PSedDiat ", 1,"x",1
"SedDiat ","fSedDiat ", 2,"x",1
"SedDiat ","VxSedDiat ", 3,"x",0
Expand Down
Binary file modified src/engines_gpl/waq/resources/process_lib/proc_def.dat
View file
Open in desktop
Binary file not shown.
Binary file modified src/engines_gpl/waq/resources/process_lib/proc_def.def
View file
Open in desktop
Binary file not shown.
233 changes: 196 additions & 37 deletions src/engines_gpl/waq/waq_process/sedim.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,9 +66,6 @@ subroutine sedim (process_space_real, fl, ipoint, increm, num_cells, &
use m_extract_waq_attribute
USE BottomSet ! Module with definition of the waterbottom segments

IMPLICIT REAL (A-H, J-Z)
IMPLICIT INTEGER (I)

REAL(kind = real_wp) :: process_space_real (*), FL (*)
INTEGER(kind = int_wp) :: IPOINT(*), INCREM(*), num_cells, NOFLUX, &
IEXPNT(4, *), IKNMRK(*), num_exchanges_u_dir, num_exchanges_v_dir, num_exchanges_z_dir, num_exchanges_bottom_dir
Expand All @@ -77,23 +74,37 @@ subroutine sedim (process_space_real, fl, ipoint, increm, num_cells, &

REAL(kind = real_wp) :: PSEDMIN

IP1 = IPOINT(1)
IP2 = IPOINT(2)
IP3 = IPOINT(3)
IP4 = IPOINT(4)
IP5 = IPOINT(5)
IP6 = IPOINT(6)
IP7 = IPOINT(7)
IP8 = IPOINT(8)
IP9 = IPOINT(9)
IP10 = IPOINT(10)
IP11 = IPOINT(11)
IP12 = IPOINT(12)
IP13 = IPOINT(13)
IP14 = IPOINT(14)
IP15 = IPOINT(15)
IP16 = IPOINT(16)
IP17 = IPOINT(17)
real(kind = real_wp), parameter :: seconds_per_day = 86400.0

integer :: ip1, ip2, ip3, ip4, ip5, ip6, ip7, ip8, ip9, ip10, &
ip11, ip12, ip13, ip14, ip15, ip16, ip17, ip18, ip19, ip20, ip21, &
in1, in2, in3, in4, in5, in6, in7, in8, in9, in10, &
in11, in12, in13, in14, in15, in16, in17, in18, in19, in20, in21
integer :: ivan, inaar, ip, ipv, ipn, ipq, iq, iseg, iflux, ikmrkv, ikmrkn, iwa1, iwa2, ikmrk1, ikmrk2, ik, iwater
real(kind = real_wp) :: conc, zersed, vsed, tau, tcrsed, delt, psed, alpha, p, pmax, maxsed, potsed
real(kind = real_wp) :: flowrate, volume, surf

IP1 = IPOINT(1) ! Concentration inorganic matter
IP2 = IPOINT(2) ! Zeroth-order flux -- at all useful?
IP3 = IPOINT(3) ! Sedimentation velocity
IP4 = IPOINT(4) ! Bottom shear stress
IP5 = IPOINT(5) ! Critical shear stress for sedimentation
IP6 = IPOINT(6) ! Depth of the segments
IP7 = IPOINT(7) ! Time step (for limiting the deposition flux)
IP8 = IPOINT(8) ! Minimum depth for sedimentation/resuspension -- obsolete!
IP9 = IPOINT(9) ! Fraction going directly to layer S2
IP10 = IPOINT(10) ! Fraction total inorganic matter (TIM) in layer S2
IP11 = IPOINT(11) ! Maximum allowable fraction (TIM) in layer S2
IP12 = IPOINT(12) ! Minimum sedimentation probability (force sedimentation) -- useful?
IP13 = IPOINT(13) ! Volume of the segments
IP14 = IPOINT(14) ! Workspace (input), used for computing maximum allowed flux
IP15 = IPOINT(15) ! (exchange) sedimentation velocity per exchange
IP16 = IPOINT(16) ! (exchange) Flow rate, used for estimating the influx
IP17 = IPOINT(17) ! (output) Sedimentation probability (used for adsorbed substances)
IP18 = IPOINT(18) ! (output) Sedimentation flux to layer S1
IP19 = IPOINT(19) ! (output) Sedimentation flux to layer S2
IP20 = IPOINT(20) ! Workspace (output), storing maximum allowed flux
IP21 = IPOINT(21) ! (additional velocity) Sedimentation velocity

IN1 = INCREM(1)
IN2 = INCREM(2)
Expand All @@ -112,15 +123,154 @@ subroutine sedim (process_space_real, fl, ipoint, increm, num_cells, &
IN15 = INCREM(15)
IN16 = INCREM(16)
IN17 = INCREM(17)
IN18 = INCREM(18)
IN19 = INCREM(19)
IN20 = INCREM(20)
IN21 = INCREM(21)

!
! Determine the mass influx to be expected for the cells near the bottom
!
! Note: this does not work properly for floating material!
!
!

do iseg = 1,num_cells
ip = ipoint(14) + (iseg - 1) * increm(14)
process_space_real(ip) = 0.0
enddo

!
! Limit the estimation to the vertical exchanges
!
DO IQ = num_exchanges_u_dir + num_exchanges_v_dir +1, num_exchanges_u_dir + num_exchanges_v_dir + num_exchanges_z_dir

IVAN = IEXPNT(1, IQ)
INAAR = IEXPNT(2, IQ)

ipv = ipoint(14) + (ivan - 1) * increm(14)
ipn = ipoint(14) + (inaar - 1) * increm(14)
ipq = ipoint(16) + (iq - 1) * increm(16)
flowrate = process_space_real(ipq)

IKMRKV = 1
IKMRKN = 1
IF (IVAN > 0 ) THEN
CALL extract_waq_attribute(1, IKNMRK(IVAN), IKMRKV)
IF ( IKMRKV /= 0 ) THEN
CALL extract_waq_attribute(2, IKNMRK(IVAN), IKMRKV)
ENDIF
ENDIF
IF (INAAR > 0 ) THEN
CALL extract_waq_attribute(1, IKNMRK(INAAR), IKMRKN)
IF ( IKMRKN /= 0 ) THEN
CALL extract_waq_attribute(2, IKNMRK(INAAR), IKMRKN)
ENDIF
ENDIF

!
! Note: including the horizontal exchanges fails, unclear why.
! Hence commented out.
!
! Handle "from" segment
!
!IF ( (IKMRKV == 0 .OR. IKMRKV == 3 ) ) THEN
! if ( inaar > 0 ) then
! ip1 = ipoint(1) + (inaar-1) * increm(1)
! else
! ! Use the "to" segment as a proxy
! ip1 = ipoint(1) + (ivan-1) * increm(1)
! endif
!
! conc = process_space_real(ip1)
!
! process_space_real(ipv) = process_space_real(ipv) + max( 0.0, -flowrate * max( 0.0, conc ) )
!
!ENDIF

!
! Handle "to" segment
!
!IF ( (IKMRKN == 0 .OR. IKMRKN == 3 ) ) THEN
! if ( ivan > 0 ) then
! ip1 = ipoint(1) + (ivan-1) * increm(1)
! else
! ip1 = ipoint(1) + (inaar-1) * increm(1)
! endif
!
! conc = process_space_real(ip1)
!
! process_space_real(ipn) = process_space_real(ipn) + max( 0.0, flowrate * max( 0.0, conc ) )
!
!ENDIF

! If the exchange is a vertical exchange, also include the settling from the cell above

if ( iq > num_exchanges_u_dir + num_exchanges_v_dir ) then
if ( ivan > 0 .and. inaar > 0 ) then
ip1 = ipoint(1) + (ivan-1) * increm(1)
ip3 = ipoint(3) + (ivan-1) * increm(3)
ip6 = ipoint(6) + (ivan-1) * increm(6)
ip13 = ipoint(13) + (ivan-1) * increm(13)
conc = process_space_real(ip1)
vsed = process_space_real(ip3) / seconds_per_day ! Conversion from m/day to m/s
depth = process_space_real(ip6)
volume = process_space_real(ip13)
if ( depth > 0.0 ) then
surf = volume / depth
else
surf = 0.0
endif
!! process_space_real(ipn) = process_space_real(ipn) + max( 0.0, -vsed * surf * max( 0.0, conc ) )
process_space_real(ipn) = max( 0.0, -vsed * surf * max( 0.0, conc ) )

ip1 = ipoint(1) + (inaar-1) * increm(1)
ip3 = ipoint(3) + (inaar-1) * increm(3)
ip6 = ipoint(6) + (inaar-1) * increm(6)
ip13 = ipoint(13) + (inaar-1) * increm(13)
conc = process_space_real(ip1)
vsed = process_space_real(ip3) / seconds_per_day ! Conversion from m/day to m/s
depth = process_space_real(ip6)
volume = process_space_real(ip13)
if ( depth > 0.0 ) then
surf = volume / depth
else
surf = 0.0
endif
!! process_space_real(ipv) = process_space_real(ipv) + max( 0.0, vsed * surf * max( 0.0, conc ) )
process_space_real(ipv) = max( 0.0, vsed * surf * max( 0.0, conc ) )

endif
endif
end do

!
! Properly scale the maximum flux
!
do iseg = 1,num_cells
ip = ipoint(14) + (iseg - 1) * increm(14)
ipv = ipoint(13) + (iseg - 1) * increm(13)
volume = process_space_real(ip13)

if ( volume > 0.0 ) then
process_space_real(ip) = seconds_per_day * process_space_real(ip) / volume
else
process_space_real(ip) = 0.0
endif
enddo

IFLUX = 0
IP1 = IPOINT(1)
IP3 = IPOINT(3)
IP6 = IPOINT(6)
IP13 = IPOINT(13)
DO ISEG = 1, num_cells

! zero output

process_space_real (IP14) = 0.0
process_space_real (IP15) = 0.0
process_space_real (IP16) = 0.0
process_space_real (IP17) = 0.0
process_space_real (IP18) = 0.0
process_space_real (IP19) = 0.0

! sedimentation towards the bottom

Expand Down Expand Up @@ -175,17 +325,18 @@ subroutine sedim (process_space_real, fl, ipoint, increm, num_cells, &
POTSED = ZERSED + (VSED * CONC) * PSED

! limit sedimentation to available mass (M/L2/DAY)
MAXSED = MIN (POTSED, CONC / DELT * DEPTH)
!!MAXSED = MIN (POTSED, CONC / DELT * DEPTH)
MAXSED = MIN (POTSED, (process_space_real(ip14) + conc / delt) * depth )

! convert sedimentation to flux in M/L3/DAY
FL(1 + IFLUX) = MAXSED * (1. - ALPHA) / DEPTH
FL(2 + IFLUX) = MAXSED * ALPHA / DEPTH
ENDIF

! Output of calculated sedimentation rate
process_space_real (IP14) = PSED
process_space_real (IP15) = MAXSED * (1. - ALPHA)
process_space_real (IP16) = MAXSED * ALPHA
process_space_real (IP17) = PSED
process_space_real (IP18) = MAXSED * (1. - ALPHA)
process_space_real (IP19) = MAXSED * ALPHA
!
ENDIF
ENDIF
Expand All @@ -206,23 +357,26 @@ subroutine sedim (process_space_real, fl, ipoint, increm, num_cells, &
IP14 = IP14 + IN14
IP15 = IP15 + IN15
IP16 = IP16 + IN16
IP17 = IP17 + IN17
IP18 = IP18 + IN18
IP19 = IP19 + IN19
end do
!
IP1 = IPOINT(1)
IP6 = IPOINT(6)
IP8 = IPOINT(8)
IP15 = IPOINT(15)
IP18 = IPOINT(18)

!.....Exchangeloop over de horizontale richting
DO IQ = 1, num_exchanges_u_dir + num_exchanges_v_dir

process_space_real(IP17) = 0.0
process_space_real(IP21) = 0.0

IP17 = IP17 + IN17
IP21 = IP21 + IN21

end do

IP13 = IP13 + (num_exchanges_u_dir + num_exchanges_v_dir) * IN13
IP15 = IP15 + (num_exchanges_u_dir + num_exchanges_v_dir) * IN15

!.....Exchangeloop over de verticale richting
DO IQ = num_exchanges_u_dir + num_exchanges_v_dir + 1, num_exchanges_u_dir + num_exchanges_v_dir + num_exchanges_z_dir + num_exchanges_bottom_dir
Expand Down Expand Up @@ -259,18 +413,18 @@ subroutine sedim (process_space_real, fl, ipoint, increm, num_cells, &
MINDEP = process_space_real(IP8 + (IVAN - 1) * IN8)
MINDE2 = process_space_real(IP8 + (INAAR - 1) * IN8)
IF (DEPTH > MINDEP .AND. DEPTH2 > MINDE2) THEN
process_space_real(IP17) = process_space_real(IP13) / 86400.
process_space_real(IP21) = process_space_real(IP15) / seconds_per_day
ELSE
process_space_real(IP17) = 0.0
process_space_real(IP21) = 0.0
ENDIF
ELSE
process_space_real(IP17) = 0.0
process_space_real(IP21) = 0.0
ENDIF

ENDIF

IP13 = IP13 + IN13
IP17 = IP17 + IN17
IP15 = IP15 + IN15
IP21 = IP21 + IN21

end do

Expand All @@ -293,6 +447,11 @@ subroutine sedim (process_space_real, fl, ipoint, increm, num_cells, &
IP15 = IPOINT(15)
IP16 = IPOINT(16)
IP17 = IPOINT(17)
IP17 = IPOINT(17)
IP18 = IPOINT(18)
IP19 = IPOINT(19)
IP20 = IPOINT(20)
IP21 = IPOINT(21)

DO IK = 1, Coll%current_size

Expand Down Expand Up @@ -336,7 +495,7 @@ subroutine sedim (process_space_real, fl, ipoint, increm, num_cells, &
ENDIF

IF (CONC > 1.E-10) THEN
process_space_real(IP17 + (IQ - 1) * IN17) = MAXSED / 86400. / CONC
process_space_real(IP21 + (IQ - 1) * IN21) = MAXSED / 86400. / CONC
ENDIF

ENDDO
Expand Down