##============================================================================== ## Advective-diffusive transport in a river ##============================================================================== tran.volume.1D <- function(C, C.up=C[1], C.down=C[length(C)], C.lat=C, F.up=NULL, F.down=NULL, F.lat=NULL, Disp, flow=0, flow.lat=NULL, AFDW = 1, V=NULL, full.check = FALSE, full.output = FALSE) { ## INPUT CHECKS N <- length(C) if (N == 0) stop("Error: C should be a vector with numeric values") if (!is.list(AFDW)) AFDW <- list(int=AFDW) if (!is.list(F.lat)) F.lat <- list(mid=rep(F.lat,length.out=N)) if (!is.list(Disp)) Disp <- list(int=rep(Disp,length.out=N+1)) if (!is.list(flow)) flow<-list(int=flow) # check later if one value/vector if (!is.list(flow.lat)) flow.lat <-list(mid=flow.lat) # check later if one value/vector if (!is.list(V)) V <- list(mid=rep(V,length.out=N)) if (!is.list(C.lat)) C.lat <- list(mid=rep(C.lat,length.out=length(C))) ## Start full check when requested if (full.check) { ## check input of grids gn <- names(AFDW) if (! "int" %in% gn) stop("error: AFDW should be a list that contains 'int', the AFDW values at the interface of the grid cells") if (is.null(AFDW$int) || length(AFDW$int) == 0) stop("error: AFDW should contain (numeric) values") if (any(AFDW$int < 0)||any(AFDW$int > 1)) stop("error: AFDW values should always range between 0 and 1") ## check input of V gn <- names(V) if (! "mid" %in% gn) stop("error: V should be a list that contains 'mid'") if (is.null(V$mid)) stop("error: the argument V should contain (numeric) values") if (any(V$mid < 0)) stop("error: teh volume V should always be larger than zero") ## check input of Disp gn <- names(Disp) if (! "int" %in% gn) stop("error: Disp should be a list that contains 'int', the dispersion coefficient at the grid cell interfaces") if (is.null(Disp$int)) stop("error: the argument Disp should be contain (numeric) values") if (any (Disp$int < 0)) stop("error: the bulk dispersion coefficient Disp should always be positive") ## check input of flow gn <- names(flow) if (! "int" %in% gn) stop("error: flow should be a list that contains 'int'") if (is.null(flow$int)) stop("error: the argument flow should contain (numeric) values") # if (any (flow$int < 0) & any (flow$int > 0)) # stop("error: the discharge flow cannot be both positive and negative within the same domain") ## check input of flow.lat gn <- names(flow.lat) if (! "mid" %in% gn) stop("error: flow.lat should be a list that contains 'mid'") # if (is.null(flow.lat$mid)) # stop("error: the argument flow.lat should contain (numeric) values") ## check input of flow gn <- names(C.lat) if (! "mid" %in% gn) stop("error: the argument C.lat should be a list that contains 'mid'") if (any (C.lat$int < 0)) stop("error: the concentration C.lat should always be positive") ## check input of flow gn <- names(F.lat) if (! "mid" %in% gn) stop("error: the argument F.lat should be a list that contains 'mid'") } # end full.check ## FUNCTION BODY: CALCULATIONS ## Calculate the discharge in each box (using the water balance) or the upstream ## flow velocity ## 1. flow.lat is NULL -> estimate flow.lat based on flow gradient (diff(flow)) if (is.null(flow.lat$mid)) { flow$int <- rep(flow$int,length.out=N+1) flow.lat$mid <- diff(flow$int) } else ## 2. flow.lat contains numeric values - check length; should be N; ## flow should contain upstream flow (one value); create flows at interface { if (!is.vector(flow.lat$mid)) flow.lat$mid <- rep(flow.lat$mid,length.out=N) if (length(flow.lat$mid) != N) stop ("flow.lat should be one number or a vector of length = N") if (!(length(flow$int) = 1)) stop ("flow should be of length = 1 if flow.lat is specified") f1 <- flow$int[1] flow$int <- c(f1,f1+cumsum(flow.lat$mid)) } ## Calculate the lateral mass input - IF OUTPUT: C is transported instead. if (is.null (F.lat$mid)) { if (is.null (C.lat$mid)) stop ("C.lat and F.lat cannot be both NULL") v1 <- flow.lat$mid F.lat$mid <- 0 if (any(flow.lat$mid>0)) { v1[v1<0] <- 0 F.lat$mid <- C.lat$mid * v1 } if (any(flow.lat$mid < 0)) { ii <- which (flow.lat$mid<0) F.lat$mid[ii] <- C[ii]*flow.lat$mid[ii] } } else { C.lat$mid <- F.lat$mid / flow.lat$mid } ## Calculate diffusive part of the mass flow F Dif.F <- as.vector(-Disp$int*diff(c(C.up,C,C.down))) ## Calculate advective part of the mass flow F ## positive flows first v1 <- flow$int Adv.F <- 0 if ( any (v1 > 0 )) { v1[v1<0] <- 0 conc <- AFDW$int *c(C.up,C) if (any (AFDW$int < 1)) conc <- conc +(1-AFDW$int)*c(C,C.down) Adv.F <- Adv.F + as.vector(v1 * conc) } else Adv.F <- 0 ## If there are negative flows: if ( any (flow$int < 0 )) { v1 <- flow$int v1[v1>0] <- 0 conc <- AFDW$int*c(C,C.down) if (any (AFDW$int < 1)) conc <- conc +(1-AFDW$int)*c(C.up,C) Adv.F <- Adv.F + as.vector(v1 * conc) } ## Assemble the total mass flow F <- as.vector(Dif.F + Adv.F) ## Impose boundary fluxes when needed if (! is.null (F.up)) F[1] <- F.up if (! is.null (F.down)) F[length(F)] <- F.down ## Calculate rate of change = Flux gradient + lateral input dC <- -diff(F)/V$mid + F.lat$mid/V$mid if (!full.output){ return (list (dC = dC, # Rate of change due to advective-diffuisve transport in each grid cell F.up = F[1], # Flux across upstream boundary ; positive = IN F.down = F[length(F)] # Flux across downstream boundary ; positive = OUT )) } else { return (list (dC = dC, # Rate of change in the centre of each grid cells flow = flow$int, flow.up = flow$int[1], flow.down = flow$int[N+1], flow.lat = flow.lat$mid, F = F, # Flux across at the interface of each grid cell F.up = F[1], # Flux across upstream boundary ; positive = IN F.down = F[length(F)], # Flux across downstream boundary ; positive = OUT F.lat = F.lat$mid)) # Flux lateral ; positive = IN } }