Which scheme for inhomogeneous convection-diffusion equation with highly variable coefficients?

I have a 1D convection-diffusion equation $sigma_t = a(x,t) sigma_{xx}+b(x,t)sigma_x+f(x,t)$ defined on the unit interval, with nonzero Neumann boundary conditions at both ends. It should be noted that the coefficients are not strictly positive, which explains the negative Peclet number below.

From what I understand, choosing an appropriate finite difference scheme depends on the Peclet number, or the ratio between the convective and diffusive coefficients, but this ratio $b(x,t)/a(x,t)$ has a strong dependence on time and varies between -700 and 1100, with the jump between these two approximate values being discontinuous in time. Here is an example plot showing this ratio at $x=0.1$, plots at other locations in space have a similar shape:

My question is, what sort of numerical scheme would be best suited to solve this problem? The region where $Pe=0$ can be discarded; the physics of the problem imply that the solution is identically zero in this region. My initial guess is to use a different scheme for each regime; i.e. when $Peapprox -700$, use one scheme, and when $Peapprox 1100$, use a different one, and use the final result of one scheme as the initial data for the next, but I’m unsure about the stability of such an approach. I’ve never seen convection-diffusion equations in the context of a negative Peclet number, so I’m unsure what scheme to use for the first regime. I’ve read elsewhere that higher-order upwinding is useful for large Peclet numbers, which might be an possible candidate for the second regime. Thoughts?