This answer is on how to implement boundary conditions (B.C.) for variable material property. My experience is that they are 2nd order accurate, but I can not find now references where such statement is proved. I hope it helps anyway, if not, it might be useful for other readers. The most closest reference is the book of Patankar on Numerical heat and mass transfer.
1.)
Let for instance $x_{i-1/2}$ is located at boundary.
If $f_i$ is located at cell center then the flux based boundary conditions are straightforward to implement, replace in your scheme the term
$$
k_{i-1/2} \frac{f_i - f_{i-1}}{\Delta x_{i-1/2}}
$$
by the term prescribed for $k \partial_x f$ in flux based boundary condition.
If you have Dirichlet B.C., then more effort is necessary. Typically you extrapolate linearly the missing (outer) value $f_{i-1}$ from given value $f_{i-1/2}$ from B.C. and from $f_i$ and then plug such $f_{i-1}$ to your scheme. In fact you might use a quadratic extrapolation depending also on $f_{i+1}$.
2.)
If $f_i$ is located in corner locations and if $x_i$ is located at boundary and you have Dirichlet B.C., you need not to apply the scheme here, the value $f_i$ is given by B.C. If you have flux-based B.C. then one typically use a "half" finite volume, so your scheme is then
$$ \frac{\frac{f_{i+1}-f_{i}}{\Delta x_{i+1/2}}k_{i+1/2}-Flux_i}{\Delta x_{i}/2}
$$
where the value $Flux_i$ is obtained from B.C.
