Only zero overland and subsurface flows to river for Irish catchment

[AquaPore]

Wflow version checks

• I have checked that this issue has not already been reported.

• I have checked that this bug exists on the latest version of Wflow.

Reproducible Example

Dear Wflow users,

I am looking to partition the river_water__volume_flow_rate between:

1. Overland flow (I guess is performed with land_surface_water__instantaneous_volume_flow_rate);
2. Subsurface flow (I guess is performed with subsurface_water__to_river_volume_flow_rate);
3. Water storage (I guess is performed with soil_layer_water_unsaturated_zone__depth);
4. Ground water (I guess is performed with subsurface_water__hydraulic_head);

The reason why I require this information is for calibration purposes and, at a later stage, to run a nutrient model.

The Wflow model runs well with [model], type=sbm. Nevertheless,

I get land_surface_water__to_river_volume_flow_rate[iX, iY, iT]=0 & subsurface_water__to_river_volume_flow_rate=0 for all cells.

I am suspecting that the problem is because I am not using model type sbm=sediment.

In the Input Toml file

[output.netcdf_grid] # ===========================================
	path = "Output_Timoleague.nc"
   compressionlevel = 0 # Compression level [0-9]

	[output.netcdf_grid.variables]
      land_surface_water__to_river_volume_flow_rate = "OverlandFlowToRiver" # [m3 / s]
      subsurface_water__to_river_volume_flow_rate = "SubsurfaceFlowToRiver" # [m3 / s]

Then I write a code to sum up the OverlandFlowToRiver and SubsurfaceFlowToRiver at every time step

WflowOutput = NCDatasets.NCDataset(Path_Wflow_NetCDF_Output)

OverlandFlowToRive = Array(WflowOutput["OverlandFlowToRiver"])
SubsurfaceFlowToRiver = Array(WflowOutput["SubsurfaceFlowToRiver"])
Precipitation    = Array(WflowOutput["Precipitation"])

Nx, Ny, Nt = size(Precipitation)

∑SubsurfaceFlowToRiver_M3 = fill(0.0, Nt)
∑OverlandFlowToRiver_M3 = fill(0.0, Nt)

for iT=1:Nt
	 for iX=1:Nx
		for iY=1:Ny
			if OverlandFlowToRiver[iX, iY, iT] !== missing
				∑OverlandFlowToRiver_M3[iT] += OverlandFlowToRiver[iX, iY, iT] * 60*60*24 # m3/s to m3/day
			end

		        if SubsurfaceFlowToRiver[iX, iY, iT] !== missing
			       ∑SubsurfaceFlowToRiver_M3[iT] += SubsurfaceFlowToRiver[iX, iY, iT] * 60*60*24 # m3/s to m3/day
			end
		end
	end 

Current behaviour

Neverthelesssum(∑OverlandFlowToRiver_M3[iT])= 0& sum(∑SubsurfaceFlowToRiver_M3[iT]) = 0

So my understanding there is a bug.

Many thanks for any help you may provide,
Joseph A.P. Pollacco

Desired behaviour

Neverthelesssum(∑OverlandFlowToRiver_M3[iT]) ≠ 0& sum(∑SubsurfaceFlowToRiver_M3[iT]) ≠ 0

Additional Context

No response

[AquaPore]

Dear @vers-w and @JoostBuitink

This post is somewhat related to Issue_895. I tried to understaand from were the issue is coming of having land_surface_water__to_river_volume_flow_rate=0.

I was searching in the code for to_river and in the function update_overland_flow_model I am questionning if the source of error is coming from the function as println(to_river) the matrix is always 0 even on the river network?

    q_av .= 0.0
    to_river .= 0.0
    qin_av .= 0.0

function update_overland_flow_model!()

"""
Update overland flow model `KinWaveOverlandFlowModel` for a single timestep `dt`. Timestepping within
`dt` is either with a fixed timestep `dt_fixed` or adaptive.
"""
function update_overland_flow_model!(
    overland_flow_model::KinWaveOverlandFlowModel,
    domain::DomainLand,
    dt::Float64,
)
    (; inwater) = overland_flow_model.boundary_conditions
    (; alpha_term, mannings_n, beta, alpha_pow, alpha) = overland_flow_model.parameters
    (; surface_flow_width, flow_length, slope) = domain.parameters
    (; q_av, qlat, qin_av, to_river) = overland_flow_model.variables
    (; adaptive) = overland_flow_model.timestepping

    @. alpha_term = pow(mannings_n / sqrt(slope), beta)
    # use fixed alpha value based flow width
    @. alpha = alpha_term * pow(surface_flow_width, alpha_pow)
    @. qlat = inwater / flow_length

    q_av .= 0.0
    to_river .= 0.0
    qin_av .= 0.0

    t = 0.0
    while t < dt
        dt_s =
            adaptive ? stable_timestep(overland_flow_model, flow_length, 0.02) :
            overland_flow_model.timestepping.dt_fixed
        dt_s = check_timestepsize(dt_s, t, dt)
        kinwave_land_update!(overland_flow_model, domain, dt_s)
        t += dt_s
    end
    q_av ./= dt
    to_river ./= dt
    qin_av ./= dt

    println(to_river)
    return nothing
end

With kindness,
Joseph A.P. Pollacco

[vers-w]

Thanks for reporting this issue.
I cannot reproduce this issue with the Wflow Moselle test model. Could you try with this model?

[AquaPore]

Dear @vers-w,

A great thanks for helping me solve this issue, where I obtained in [output.netcdf_grid] that land_surface_water__to_river_volume_flow_rate[iX, iY, iT] =0 for IrishCatchment.

Thanks for pointing me to the rtight direction and I would like to confirm that there is no issue with Wflow Moselle . I modified the input toml such that it has the same parameters as the IrishCatchment such that in this test it models impermeable layers Issue_895 therfore the log.txt of Ireland and Moselle are quasi identical except that cell_length_in_meter__flag = true for IrishCatchment and false for Moselle.

My investigation

It is beyond my understanding why for IrishCatchment and Moselle I obtain sum(to_river ) >> 0 which I implemented in the following section of the module src/routing/surface/surface_kinwave.jl but land_surface_water__to_river_volume_flow_rate[iX, iY, iT] =0 .

function kinwave_land_update!()

surface_kinwave.jl

function update_overland_flow_model!()

update_overland_flow_model!

function update_lateral_inflow!(

update_lateral_inflow
println(sum(overland_flow.variables.to_river))

What tests will you suggest I perform?

What other functions does overland_flow.variables.to_river pass through that I may want to investigate

Toml model

[model] # ===========================================
type = "sbm" # <sbm> || <sbm_gwf> || <sediment>
cell_length_in_meter__flag = true # <false> Unit cell length of raster cells in x and y direction in lat/lon degree; <true> in meter
cold_start__flag = false # Cold state does not include values of storage at time=0
glacier__flag = false # Include glacier modelling
reservoir__flag = false # Include reservoir modelling
snow__flag = false # Include snow modelling
snow_gravitational_transport__flag = false # include lateral snow transport in the model, default is false
soil_infiltration_reduction__flag = false # <soil_infiltration_reduction> reduction of conductivity due to ice
water_mass_balance__flag = false # Compute water balance error
kinematic_wave__adaptive_time_step_flag = false # Enable kinematic wave adaptive (internal) timestepping
land_streamorder__min_count = 5 # Multi-Threading: Minimum stream order to delineate subbasins for river domain, default is 6 (for multi-threading computing purposes)
river_streamorder__min_count = 6 # Multi-Threading: Minimum stream order to delineate subbasins for land domain, default is 5 (for multi-threading computing purposes)
soil_layer__thickness = [100, 300,800 ] 
saturated_hydraulic_conductivity_profile = "exponential" # <kh_profile_type>; [char]  |"layered"|,

Many thanks for any help you may provide,
With kindness,
Joseph A.P. Pollacco

[vers-w]

Maybe you can share your complete model setup (including forcing) so we can investigate the issue?

[AquaPore]

Maybe you can share your complete model setup (including forcing) so we can investigate the issue?

Thanks @vers-w for your great kindness and you amazing support of DELTARES, the files of the small Irish catchment had been send to your email.
With kindness,
Joseph A.P. Pollacco

[vers-w]

Thanks for providing the model files. I did investigate the issue, and the to_river variables are computed correctly internally, however there is an issue with writing these variables as output.

[AquaPore]

Dear @vers-w,

I want to thanks DELTARES for your amazing support.

I suspect that the error is caused by NCDatasets.jl how it deals with NaN data translated into missing . NCDatasets has new experimental functions on how it handles missing data, which might be useful
to upgrade from "0.13.2, 0.14" to the latest version is v0.14.15.

With kindness,
Joseph A.P. Pollacco

[vers-w]

I suspect that the error is caused by NCDatasets.jl how it deals with NaN data translated into missing .

Thanks for the suggestion, it is not caused by NCDatasets. The active indices for the to_river variables (to map the model vector to the 2D model domain) are not correct, this needs to be fixed in v1.0 as well as in the Master branch.