I'm trying to understand an approximate implicit solution to a nonlinear diffusion equation from Perron (2011). The equation in question is
$$\frac{\partial z(x,t)}{\partial t} = \frac{\partial}{\partial x}\left(\dfrac{k}{1-\left(\dfrac{1}{S_c}\dfrac{\partial z(x,t)}{\partial x}\right)^m} \cdot\frac{\partial z(x,t)}{\partial x}\right) \tag{1} $$
where $z$ is the time $t$ and space $x$ dependent elevation. $S_c$, $m$, $k$ are some constants. $z(x,t)$ is abbreviated as $z_i^n$ where $i$ denotes a discrete point in space and $n$ a point in time. Somewhat confusingly, the approximations for the derivatives of $z$ with respect to $x$ are written as $z_x$ and $z_{xx}$: $$ \frac{\partial z_i^n}{\partial x} \approx z_x = \frac{z_{i+1}^n - z_{i-1}^n}{2\Delta x} $$ $$ \frac{\partial^2 z_i^n}{\partial x^2} \approx z_{xx} = \frac{z_{i+1}^n-2z_i^n+z_{i-1}^n}{(\Delta x) ^2} $$
Thus, we can write an approximation of the diffusion equation:
$$\frac{\partial z(x,t)}{\partial t} \approx k\left[\dfrac{z_{xx}}{1-\left(S_c^{-1}z_{x}\right)^m} + \dfrac{mz_x^m z_{xx}}{S_c^m\left(1-\left(S_c^{-1}z_x\right)^m\right)^2}\right] \label{nonlinDiff}\tag{2} $$
Now, Perron (2011) proposes and approximate implicit solution and states "I evaluate the flux divergence analytically, and then perform the Taylor series expansion on the finite difference approximation for the flux divergence". From my understanding, the term in big brackets in equation 1 is the flux term, and replacing the $z_{xx}$ and $z_x$ in equation 2 with the actual partial derivatives would be the analytical flux divergence. Thus, equation 3 is a finite difference approximation of the flux divergence.
What I don't understand: How do you calculate the (first order) Taylor series expansion of the right hand term in equation 2? Perron formulates the following equation:
$$ \frac{z_i^{n+1}-z_i^n}{\Delta t} = f(z_i^n) + \sum_m\left[\frac{\partial f(z_i^n)}{\partial z_m^n}(z_m^{n+1}-z_m^n)\right] \tag{3} $$
where $m = i-1,i,i+1$. $f(z_i^n)$ is defined as the right hand side of equation 2. Perron then introduces another notation:
$$ F_k = \frac{\partial f(z_i^n)}{\partial z_k^n}. $$
And actually gives the terms for $F_{i-1}$, $F_i$ and $F_{i+1}$ in the appendix:
$$ F_i = -2ac\frac{1+2bcz_x^2}{(\Delta x)^2} $$
$$ F_{i+1} = ac\frac{1+2bcz_x^2}{(\Delta x)^2} + abc^2\frac{(1+2bcz_x^2)z_xz_{xx}}{\Delta x} + 2abc^2\frac{(1+bcz_x^2)z_xz_{xx}}{\Delta x} $$
$$ F_{i-1} = ac\frac{1+2bcz_x^2}{(\Delta x)^2} - abc^2\frac{(1+2bcz_x^2)z_xz_{xx}}{\Delta x} - 2abc^2\frac{(1+bcz_x^2)z_xz_{xx}}{\Delta x} $$
where $a = k$, $b = S_c^{-2}$, $c = \left(1-bz_x^2\right)^{-1}$. Note that this is the solution to the case $m = 2$. I was able to implement this scheme in Python and it works, but I want to generalize it to any $m$. That's why I'm trying to understand how these $F_k$ terms are calculated. Maybe I got lost in all the notational conventions or I just don't understand partial derivatives well enough, but what does a term like
$$ \frac{\partial f(z_i^n)}{z_{i-1}^n}~\mathrm{or}~\frac{\partial f(z_i^n)}{z_{i}^n} $$
even mean? In my understanding, $f$ depends on $z$, which in turn depends on $t$ and $x$. So is Perron taking the partial derivative of $f$ with respect to another function? How would you calculate that?