How should the FV scheme be implemented?

[HannoSpreeuw]

According to Boudreau's formula 8.72 the central differenced first order spatial derivative should be replaced by his formula, following Fiadeiro & Veronis (FV) "On weighted-mean schemes for the finite-difference approximation to the advection-diffusion equation" (1977).

For the first term on the rhs of equations 42 from L'Heureux it is clear how the first order spatial derivative should be applied, i.e. following Boudreau. But I am not sure how that works with the second term on the rhs of equations 42 from L'Heureux. If we apply the chain rule to the second term of equations 42 we again get first order spatial derivatives, i.e. a product of the spatial derivative of the porosity times the spatial derivative of the concentration of one of the ions (plus a term with a laplacian, similar to the first term of 8.72). Should the FV way of calculating first order spatial derivatives be applied to both factors in that product? That is the way this has been done in marlpde, but not in rhythmite. It feels weird to calculate spatial derivatives in two different ways in the same equation.

[csummers25]

In formula in Boudreau, the FV scheme is only applied to the advection term, the diffusion term is just a second order centred difference. I think it makes sense therefore to only apply the FV scheme to the advective derivatives, not those from the diffusion term, even though the diffusion term in equation 42 has a different form, since the idea of the FV scheme is to add numerical diffusion to the advective terms to improve their stability.

[EmiliaJarochowska]

Hanno, I am not 100% sure what to answer, other than my understanding aligns with that @csummers25 wrote: the scheme's purpose is to add stability to advection when the Peclet number varies a lot. And Fortran, for instance, uses blending between upwind and Crank-Nicolson, rather than central differencing, for the diffusion terms, so that is already a difference between the implementations, and of course Fortran implements the F-V scheme incorrectly, it using it three different sigmas.

That's not an answer to your question, other than I am in favor of limiting the use if this scheme and the argument for using similar approaches between implementations doesn't hold anyway without either using C-N in rhytmite or central differencing in Fortran, neither of which seems needed, since rhytmite is stable even with central differencing. So I cannot give you a satisfactory answer, but based on @csummers25 email, we can do without F-V in rhytmite runs so there will be no agreement between marlpde and rhytmite with that respect anyways. I tried to calculate if the diffusion terms will not end up being the same if F-V is used but I'm not confident using just pen and paper for that.

Thinking about this I recalled this old paper that totally roasted F-V and for that is worth browsing through: https://www.osti.gov/servlets/purl/5603843-1fIMeI/ (see page 7 of the pdf)
Gresho, P. M., & Lee, R. L. (1981). Don't suppress the wiggles—they're telling you something!. Computers & Fluids, 9(2), 223-253.

[HannoSpreeuw]

Thanks @EmiliaJarochowska for that link to the paper by Gresho and Lee and their remark about F-V, i.e. that it can fail miserably for time-dependent problems.

I tried to find out if F-V still holds for anything more complex than a default advection-diffusion equation.

This is the case for equations 42, these are slightly more complex than default, because you get an extra term when you apply the derivative in the second term of the rhs, because of the chain rule.

But perhaps 43 is also slightly more complex than a default advection-diffusion equation, because you have $-\frac{\partial (W \cdot \phi)}{\partial x}$ and not $-W \cdot \frac{\partial \phi}{\partial x}$. Also, wrt to the second term on the rhs of 43, $d_{\phi}$ is depends on time and depth.

Fiadeiro and Veronis (1977) write (paragraph 3, page 516):

In more complex situations, ones with variable coefficients and/or more than one dimension, 
there is no general technique for obtaining an exact equivalent equation.

So perhaps the F-V scheme is not applicable to equations 42 and 43.

Something similar is written in paragraph 4 of "Finite difference approximations to the advection-diffusion equation" by Daniel G. Wright:

The numerical schemes discussed in Sections 2 and 3 give exact results for the solutions to (1) and (11) 
no matter how coarse the resolution is chosen to be. However, they are clearly not exact if u and K 
are spatially varying.