I am trying to implement a 1D k-$\varepsilon$ turbulent model for a turbidity current, hence the conservation equation for $c$. I'm solving for the variables $u,c,k$ and $\varepsilon$. The remaining terms are parameters.$$\newcommand{\cmu}{c_{\mu}}\newcommand{\cmuc}{c_{\mu c}}\newcommand{\dervz}[1]{\frac{\partial #1}{\partial z}}\newcommand{\eddyvis}{\cmu \frac{k^2}{\varepsilon}}\newcommand{\eddyvisC}{\cmuc \frac{\hat{k}^2}{\hat{\varepsilon}}}$$$$\frac{\partial u }{\partial t} = \frac{\partial}{\partial z} \left[ \left( \nu + \cmu \frac{k^2}{\varepsilon} \right) \frac{\partial u }{\partial z} \right]+ \frac{\rho_s - \rho_w}{\rho_w} c g_x$$$$ \frac{\partial c}{\partial t} = v_s \dervz{c} + \frac{\partial}{\partial z} \left[ \left( \nu + \cmuc \frac{k^2}{\varepsilon} \right) \dervz{c} \right] $$$$ \frac{\partial k}{\partial t} =  \frac{\partial}{\partial z}  \left[ \left( \nu + \frac{\cmu}{\sigma_k} \frac{k^2}{\varepsilon} \right) \dervz{k} \right] + \eddyvis \left(\dervz{U}\right)^2 - \varepsilon   + \frac{\rho - \rho_w}{\rho_w} g_z \frac{1}{\sigma} \eddyvisC \dervz{c}$$$$ \frac{\partial \varepsilon}{\partial t} =  \frac{\partial}{\partial z}  \left[ \left( \nu + \frac{\cmu}{\sigma_{\varepsilon}} \frac{k^2}{\varepsilon} \right) \dervz{\varepsilon} \right]   + c_{\varepsilon 1} \frac{\varepsilon}{k}   \left( \eddyvis \left(\dervz{u}\right)^2  + c_{\varepsilon 3}   \frac{\rho - \rho_w}{\rho_w} g_z \frac{1}{\sigma} \eddyvisC \dervz{c} \right) - C_{\varepsilon 2} \frac{ \varepsilon^2}{k}$$So far, I am using finite difference. Upwind scheme for the first order derivatives. Because the paramter $v_s$ is positive, concentration travels downwards at a speed $v_s$, I'm implementing the upwind scheme as:$$\frac{c^{i+1}_n - c^{i}_n}{\Delta z}$$for the node $i$. I also use the same upwind scheme for all the first-order derivative terms, i.e. also the terms inside the diffusion terms, which are nonlinear. For the linear diffusion terms, I apply second order central difference. For the time integration, I use a Crank-Nicholson method treating all terms implicitly. The boundary conditions are complex. I'm using a standard wall function, using the log law at a node far from the wall. Using the log law, I also impose a gradient boundary condition because the log law adds the shear velocity as a unknown, hence the additional equation/condition. For $k$ and $\varepsilon$ I use conditions that depend on the shear velocity. That's for one side of the domain, at the other side, sufficiently far, I use zero value conditions.I'm not seeing much success with this implementation, should I use different discretization schemes for the diffusion and the advective terms? Should I use a different method altogether such as Finite Volume? What confuses me are the diffusion terms. Which discretization should I use there?With regards to the time discretization. Should I treat some terms explicitely? The source terms are nonlinear as well.Thanks in advance.EDIT: Response to the first commentSorry I was not clear there. The model blows up for certain initial values of the concentration (too much acceleration), although that might have to do with the wall function implementation I think. Another issue that I had  were the initial conditions for $k$ and $\varepsilon$. I'll explain. The initial conditions for the velocity and concentration do not span the entire domain. They have positive values only for one tenth of the lower part of the domain, where I implement the wall function boundary condition, eventually, they will diffuse upwards. However, they will not do so if the values for $k$ and $\varepsilon$ are zero above where $u$ and $c$ are. This is because they are in the diffusion term and I think because of the upwind scheme, which assumes the wave travels downwards.
1D k-epsilon turbulence model in a turbidity current
finite difference;advection diffusion