For the past few months, I've been working on a similar problem. I only solve the steady-state potential equation for the moment, assuming homogeneous, steady-state ion concentrations.
Scaling your equations is going to be a rudimentary preconditioner. I found that I got good results using algebraic multigrid (AMG) for preconditioning instead, since my problem reduces to a coupled reaction-diffusion equation (i.e., a Poisson problem with nonlinear source term). You might try incomplete LU (ILU) to start, because it is a robust, slower preconditioner, and if that works, then try other preconditioners like AMG.
Using a direct solver might be helpful, but it really depends on the size of your problem. Even sparse LU decompositions will exhaust available memory due to fill-in for problems that are sufficiently large (that point occurs somewhere in the neighborhood of 100,000-1,000,000 variables). If your Jacobian matrix is rank-deficient (as was the case in my problem, due to a nonconductive filler material in some cells), then direct methods will also have problems.
A single higher quality guess for the nonlinear solver won't help. Your problem is transient, so you'll need a guess for every time step, if you're using an implicit method (and given that it's a diffusion problem, you probably should be, unless you want your time steps to be severely limited by the diffusion stability limit). Established libraries that implement numerical integration methods (i.e., time-steppers) do a good job of constructing initial guesses for each time step, and if the guesses are not great (they're typically derived from a problem-independent heuristic), the time steps will be reduced to preserve accuracy.
Quad precision might be helpful, if your problem defies preconditioning. There aren't many libraries I'm aware of that use quad precision. PETSc does, but I believe if you configure with quad precision, you're limited to PETSc-native solvers and preconditioners.
It's possible. You should test your assembled Jacobian against a finite difference approximation for debugging purposes. If you are using a finite difference approximation in your actual simulation, you will roughly halve the number of significant figures in your problem, which could severely erode your accuracy because your problem is ill-conditioned; that said, it doesn't seem like you're using that approximation in your simulation.