Julia ModelingToolkit Method for Dx(var) * Dx(var)

I am trying to use the Julia language to solve a system of partial differential equations. I found the modeling tool kit to be quite helpful. I can solve the unsteady heat equation for multiple species and for Neumann and Dirichlet boundary conditions, see MWE.

However, I am looking to solve more complicated problems, such as problems with electromigration, where the flux is

N_i = -D_i\nabla(c_i) - z_i u_i F c_i \nabla(\Phi_2)

From the MWE, I think I understand how to solve a basic PDE problem

MWE:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets
using Plots

# parameters, variables, derivatives
@parameters t x
@variables c1(..) c2(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

# Define flux
function diffusive_flux_func(diff_val, u)
    # Example
    # diffusive_flux_1 = diffusive_flux_func(diff_1, c(t,x))
    return -diff_val * Dx(u)
end

diff_1 = 1.0 / 10.0
diff_2 = 2.0 / 10.0

diffusive_flux_1 = diffusive_flux_func(diff_1, c1(t, x))
diffusive_flux_2 = diffusive_flux_func(diff_2, c2(t, x))

source_term_1 = 0.0
source_term_2 = 0.0

conservation_eq_1 = Dt(c1(t, x)) ~ -Dx(diffusive_flux_1) + source_term_1
conservation_eq_2 = Dt(c2(t, x)) ~ -Dx(diffusive_flux_2) + source_term_2

# governing equations
eq = [conservation_eq_1, conservation_eq_2]

# fluxes at x = 0.0
flux_c1_x0 = diffusive_flux_func(diff_1, c1(t, 0))
flux_c2_x0 = diffusive_flux_func(diff_2, c2(t, 0))

# fluxes at x = xmax
flux_c1_xmax = diffusive_flux_func(diff_1, c1(t, 1.0))
flux_c2_xmax = diffusive_flux_func(diff_2, c2(t, 1.0))

bcs = [
    c1(0, x) ~ 0.0,  # initial condition
    flux_c1_x0 ~ 0.0,  # bc at x = 0 - Neumann
    c1(t, 1) ~ 1.0,   # bc at x = 1 - Dirichlet

    c2(0, x) ~ 1.0,
    c2(t, 0) ~ 0.0,
    flux_c2_xmax ~ 0.0,
]

# space and time domains
domains = [
    t ∈ Interval(0.0, 1.0),
    x ∈ Interval(0.0, 1.0)
]

# PDE System
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [c1(t, x), c2(t, x)])

# Discretization
dx = 0.02
order = 2
discretization = MOLFiniteDifference([x => dx], t)

# Convert PDE into an ODE
prob = discretize(pdesys, discretization)

# solve the ODE problem
sol = solve(prob, Tsit5(), saveat=0.1)
display(sol)

# Plot results
# Magic indexing with symbolic variables (custom getindex methods through julia package)
x_vals = sol[x]
t_vals = sol[t]
c1_vals = sol[c1(t,x)]
c2_vals = sol[c2(t,x)]

plt = plot(layout=(1,2), size=(1000, 400))

for i in eachindex(t_vals)
    plot!(x_vals, c1_vals[i, :], label="t=$(t_vals[i])", subplot=1)
    plot!(x_vals, c2_vals[i, :], label="t=$(t_vals[i])", subplot=2)
end

display(plt)
# save plot
savefig(plt, "diffusion_2species.png")

What I tried:

# Constants
Rigc = 8.314
Fconst = 96485
Temp = 298

# parameters, variables, derivatives
@parameters t x
@variables c_cat(..) c_an(..) Φ₂(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

# Define flux
# diffusion migration flux
function diffusion_migration_flux_func(D_i, z_i, c_i, Φ₂)
    # Example:
    # flux_cat = diffusion_migration_flux_func(diff_cat, z_cat, c_cat(t, x), Φ₂(t, x))
    u_i = D_i / (Rigc * Temp)   # Nernst-Einstein relation

    # Nᵢ = -Dᵢ∇(cᵢ) - zᵢ⋅uᵢ⋅F⋅(cᵢ∇Φ₂)
    flux_i = -D_i * Dx(c_i) - z_i * u_i * Fconst * c_i * Dx(Φ₂)
    return flux_i
end

diff_cat = 1.0e-3
z_cat = +1.0
u_cat = diff_cat / (Rigc * Temp)

diff_an = 2.0e-3
z_an = -1.0
u_an = diff_an / (Rigc * Temp)

flux_cat = diffusion_migration_flux_func(diff_cat, z_cat, c_cat(t, x), Φ₂(t, x))
flux_an = diffusion_migration_flux_func(diff_an, z_an, c_an(t, x), Φ₂(t, x))

source_term_cat = 0.0
source_term_an = 0.0

# how I want to write the equations: control-volume formulation
conservation_eq_cat = Dt(c_cat(t, x)) ~ -Dx(flux_cat) + source_term_cat
conservation_eq_an = Dt(c_an(t, x)) ~ -Dx(flux_an) + source_term_an

# differential form
# conservation_eq_cat = Dt(c_cat(t, x)) ~ (diff_cat * Dxx(c_cat(t,x)) 
#                                             + z_cat * u_cat * Fconst * Dx(c_cat(t,x)) * Dx(Φ₂(t,x))
#                                             + z_cat * u_cat * Fconst * c_cat(t,x) * Dxx(Φ₂(t,x))
#                                             )
# conservation_eq_an = Dt(c_an(t, x)) ~ (diff_an * Dxx(c_an(t,x)) 
#                                             + z_an * u_an * Fconst * Dx(c_an(t,x)) * Dx(Φ₂(t,x))
#                                             + z_an * u_an * Fconst * c_an(t,x) * Dxx(Φ₂(t,x))
#                                             )
electrolyte_potential_eq = 0.0 ~ Dxx(Φ₂(t, x))

# governing equations
eqs = [
    conservation_eq_cat, 
    conservation_eq_an, 
    electrolyte_potential_eq,
    ]

# fluxes at x = 0.0
flux_c_cat_x0 = diffusion_migration_flux_func(diff_cat, z_cat, c_cat(t, 0.0), Φ₂(t, 0.0))
flux_c_an_x0 = diffusion_migration_flux_func(diff_an, z_an, c_an(t, 0.0), Φ₂(t, 0.0))

# fluxes at x = xmax
flux_c_cat_xmax = diffusion_migration_flux_func(diff_cat, z_cat, c_cat(t, 1.0), Φ₂(t, 1.0))
flux_c_an_xmax = diffusion_migration_flux_func(diff_an, z_an, c_an(t, 1.0), Φ₂(t, 1.0))

bcs = [
    c_cat(0, x) ~ 1.0,  # initial condition
    c_cat(t, 0.0) ~ 0.9,  # bc at x = 0 - Neumann
    c_cat(t, 1) ~ 1.0,   # bc at x = 1 - Dirichlet

    c_an(0, x) ~ 1.0,
    c_an(t, 0) ~ 0.9,
    c_an(t, 1.0) ~ 1.0,

    Φ₂(0, x) ~ 0.0,
    Φ₂(t, 0.0) ~ 0.0,
    Φ₂(t, 1.0) ~ 0.0,
]

# space and time domains
domains = [
    t ∈ Interval(0.0, 0.01),
    x ∈ Interval(0.0, 1.0)
]

# PDE System
@named pdesys = PDESystem(eqs, bcs, domains, [t, x], [c_cat(t, x), c_an(t, x), Φ₂(t, x)])

# Discretization
dx = 0.02
order = 2
discretization = MOLFiniteDifference([x => dx], t)

# Convert PDE into an ODE
prob = discretize(pdesys, discretization)

# solve the ODE problem
# Tsit5 not able to use mass matrices
# FBDF() can be used for DAE systems
sol = solve(prob, FBDF(), saveat=0.001)
display(sol)

# Plot results
# Magic indexing with symbolic variables (custom getindex methods through julia package)
x_vals = sol[x]
t_vals = sol[t]
c_cat_vals = sol[c_cat(t,x)]
c_an_vals = sol[c_an(t,x)]
Φ_2_vals = sol[Φ₂(t,x)]

plt = plot(layout=(1,3), size=(1400, 400))

for i in eachindex(t_vals)
    plot!(x_vals, c_cat_vals[i, :], label="t=$(t_vals[i])", subplot=1)
    plot!(x_vals, c_an_vals[i, :], label="t=$(t_vals[i])", subplot=2)
    plot!(x_vals, Φ_2_vals[i, :], label="t=$(t_vals[i])", subplot=3)
end

display(plt)
# save plot
savefig(plt, "Electromigration_flux.png")

The solution is not physical and for some reason additional dependent variables are added:

Dependent variables:
    Φ₂(t, x): (11, 51) sized solution
    c_cat(t, x): (11, 51) sized solution
    var"⟦-0.002Differential(x)(c_an(t, x)) + 0.07788673749945511Differential(x)(Φ₂(t, x))*c_an(t, x)⟧"(t, x): (11, 51) sized solution
    var"⟦-0.001Differential(x)(c_cat(t, x)) - 0.03894336874972756Differential(x)(Φ₂(t, x))*c_cat(t, x)⟧"(t, x): (11, 51) sized solution
    c_an(t, x): (11, 51) sized solution

What is a method for solving a PDE using Julia's ModelingToolkit when there are terms such as Dx(c_cat(t, x)) * Dx(Phi_2(t, x)). This seems to be the crux of the problem for me.