First some notation
$$\nabla \times \vec B \rightarrow \epsilon_{ijk}\nabla_j B_k$$
$$\nabla \cdot \vec B \rightarrow \nabla_i B_i$$
$$\nabla B \rightarrow \nabla_i B$$
Now, to your problem,
$$\nabla \cdot(\nabla \times \vec V)$$
writing it in index notation
$$\nabla_i (\epsilon_{ijk}\nabla_j V_k)$$
Now, simply compute it, (remember the Levi-Civita is a constant)
$$\epsilon_{ijk} \nabla_i \nabla_j V_k$$
Here we have an interesting thing, the Levi-Civita is completely anti-symmetric on i and j and have another term $\nabla_i \nabla_j$ which is completely symmetric: it turns out to be zero.
$$\epsilon_{ijk} \nabla_i \nabla_j V_k = 0$$
Lets make the last step more clear. We can always say that $a = \frac{a+a}{2}$, so we have
$$\epsilon_{ijk} \nabla_i \nabla_j V_k = \frac{1}{2} \left[ \epsilon_{ijk} \nabla_i \nabla_j V_k + \epsilon_{ijk} \nabla_i \nabla_j V_k \right]$$
Now lets interchange in the second Levi-Civita the index $\epsilon_{ijk} = - \epsilon_{jik}$, so that
$$\epsilon_{ijk} \nabla_i \nabla_j V_k = \frac{1}{2} \left[ \epsilon_{ijk} \nabla_i \nabla_j V_k - \epsilon_{jik} \nabla_i \nabla_j V_k \right]$$
Now we can just rename the index $\epsilon_{jik} \nabla_i \nabla_j V_k = \epsilon_{ijk} \nabla_j \nabla_i V_k$ (no interchange was done here, just renamed).
$$\epsilon_{ijk} \nabla_i \nabla_j V_k = \frac{1}{2} \left[ \epsilon_{ijk} \nabla_i \nabla_j V_k - \epsilon_{ijk} \nabla_j \nabla_i V_k \right]$$
We can than put the Levi-Civita at evidency,
$$\epsilon_{ijk} \nabla_i \nabla_j V_k = \frac{\epsilon_{ijk}}{2} \left[ \nabla_i \nabla_j V_k - \nabla_j \nabla_i V_k \right]$$
And, because V_k is a good field, there must be no problem to interchange the derivatives $\nabla_j \nabla_i V_k = \nabla_i \nabla_j V_k$
$$\epsilon_{ijk} \nabla_i \nabla_j V_k = \frac{\epsilon_{ijk}}{2} \left[ \nabla_i \nabla_j V_k - \nabla_i \nabla_j V_k \right]$$
And, as you can see, what is between the parentheses is simply zero.
The problem is that you used the same indices to sum over the elements of the first and second $A\times B$ in your dot product. The product should actually read
$$ \mid \; \vec{A} \times \vec{B} \mid^2 \; = (\vec{A} \times \vec{B})\cdot(\vec{A} \times \vec{B}) = \epsilon_{ijk}a_j b_k \hat{e}_i \cdot\; \epsilon_{pqr}a_q b_r \hat{e}_p$$
From here on you can use the orthogonality of the $\hat{e}_i$s and then use the following relation to simplify your problem
\begin{align}\varepsilon _{ijk}\varepsilon _{imn}=\delta _{jm}\delta _{kn}-\delta _{jn}\delta _{km}\end{align}
Best Answer
I asked in comments if an answer involving abstract index notation would be welcome, and apparently it was.
Manifolds
Because you're talking about $\nabla$, let's jump into the deep end with manifolds. You have some space of points $\mathcal M$ and some field of numbers $\mathbb F$ (usually the real numbers $\mathbb R$ or the complex numbers $\mathbb C$); it becomes a "manifold" as we start to endow it with a set of "allowed" scalar fields $\mathcal S \subseteq (\mathcal M \to \mathbb F)$ which we think of as "smooth functions." To make that point more precise, one can imagine the existing smooth functions among the numbers which we understand much better, from $\mathbb F^n \to \mathbb F$. This space $C^\infty(\mathbb F^n,~\mathbb F)$ of these functions can be used as a closure relation, with the following category diagram,
Note that this is an abstract diagram and there are circles implied for $n=0,1,2,\dots$. So on the left we have these $n$-tuples of scalar fields $\mathcal S^n$; these would normally go from $\mathcal M^n \to \mathbb F^n$ but we are applying all of the scalar fields to the same source point, so there is an implied function $\mathcal M\to\mathcal M^n$ in there as well. On the right we have these smooth functions, the "category diagram" says that composing any $n$-tuple of scalar functions with any of these smooth functions of the appropriate $n$, yields another scalar field, another of these top-arrows. And the point is, every smooth function from $\mathbb F^n$ to $\mathbb F$ now enjoys a sort of dual-status, on the one hand it is a conventional smooth function; on the other hand it can be "lifted" into a function from $\mathcal S^n \to \mathcal S,$ and we say that these scalar fields are smooth precisely because they're closed under all of these lifted functions. For short, let me call such a function an $n$-functor and denote it with square braces when applied to scalar fields or parentheses when applied to numbers, so that $$f[s_1,~s_2,~s_3](p) = f\big(s_1(p),~s_2(p),~s_3(p)\big).$$
I want to call out a couple of these functions in particular: addition of scalar fields is now defined because $(+)$ is a 2-functor, and pointwise multiplication of scalar fields is now defined because $(\cdot)$ is a 2-functor. We also have scalar multiplications and scalar additions by any constants in $\mathbb F$ because those are 1-functors, and we have constant fields because those are 0-functors. I won't use square braces for any of those. Finally it turns out that one can derive a natural topology from the set of fields you use (this is one of the reasons you don't want to just use $\mathcal S = (\mathcal M \to \mathbb F),$ you get the discrete topology on the space) -- for this it's helpful that the bump function is a 1-functor.
What makes this thing a $D$-dimensional manifold is the axiom that for every point $p\in\mathcal M$ there is some open set containing that point where $D$ of these scalar fields, "local coordinate fields", can be used to (a) distinguish points within that open set, and (b) represent any other scalar field as a $D$-functor applied to these coordinate fields. There are some other necessary axioms too, like if a scalar field is piecewise defined on a patchwork of all of these spaces consistently then it should also be in $\mathcal S$, but let's skip those for brevity.
Vector fields
There now exists a nice definition of the space of vector fields on $\mathcal M$ but it is a bit abstract: it is the space of directional derivatives, which are called derivations, on the manifold. Formally: $\mathcal V$ is the subset of functions from $\mathcal S \to \mathcal S,$ which obey the Leibniz law: if $V$ is in this subset then its action on any $n$-functor is given by $$ V\big(f[s_1,\dots s_n]\big) = \sum_{i=1}^n f_{(i)}[s_1,\dots s_n]\cdot V(s_i), $$ where $f_{(i)}$ is the derivative of $f$ with respect to its $i^\text{th}$ argument. (There's a lot of machinery going on here! I'm saying: go back to what $f$ is as a function from numbers to numbers, take the derivative, the resulting function is an $n$-functor, apply it to these scalar fields as such, then multiply these pointwise with $V$ applied to the individual scalar fields, then sum everything up.)
Why would this be a "vector field"? Well, go back to our coordinate axiom: every scalar field is a $D$-functor of the coordinate fields on the open set. This means that on this open set, given the $D$ "components" $v_i = V(c_i)$ (which are also scalar fields), the operation of $V$ on a scalar field $s$ is uniquely given by $V s = \sum_i v_i\cdot s_{(i)}.$ So those components entirely define the vector, and in fact I think you can use (by the patchwork axiom) any $D$-tuple of scalar fields to create one of these.
The space of vector fields $\mathcal V$ is, perversely, not quite a vector space over the field $\mathcal S$, and this is because $\mathcal S$ violates the field axioms: you can have one function which is zero on the top half of the sphere, and one function which is zero on the bottom half of the sphere, and multiply them to get the $0$-element that is zero on the whole sphere: and the field axioms prohibit divisors of zero. (I guess more directly: each of those functions has no multiplicative inverse, but it is not the zero element.) Instead we have to say that $\mathcal V$ is a module over the commutative ring $\mathcal S.$
Tensor algebra
Now that we have the vector fields $\mathcal V$ we invent the covector space $\bar {\mathcal V}$, which is the space of linear maps from vectors to scalar fields, $\bar{\mathcal V} = \operatorname{Hom}(\mathcal V, \mathcal S).$ This space also is a module over $\mathcal S,$ with scalar multiplication meaning "multiply the scalar-field output of this covector, pointwise, by the given scalar field" and addition meaning "add the scalar-field outputs of these two covectors." In fact for all pairs of natural numbers, there is a module of multilinear operators from $m$ covectors and $n$ vectors to scalar fields, the space of $[m, n]$-tensor fields $$\mathcal T[m, n] = \operatorname{Hom}(\mathcal V^n\times\bar{\mathcal V}^m, \mathcal S).$$ For dealing with these, there are many notations which have been invented around the world. I'm now going to show you the one known as abstract index notation.
Abstract indices
The idea is that we create copies of the $[m, n]$ tensor space for any set of $m+n$ distinct symbols and denote it accordingly, $\mathcal T^{abc}_{de}$ being one copy of the space of $[3, 2]$-tensors. Every element of this space must also be labeled with the appropriate symbols, and for individual elements they may also need be in some particular order (for the tensor-space, the top and bottom symbols are order-independent).
This reveals for us two things: a family of outer products, for example one of these outer products maps $\mathcal T^{ab}_e \times \mathcal T^c_d \mapsto \mathcal T^{abc}_{de}$. The meaning of this is pretty straightforward, I think: the multi-linear map $A^{ab}_e ~ B^c_d$ takes as argument the two vectors $u^d, ~ v^e$ and the three covectors $r_a, s_b, t_c$ and produces the two scalar fields $A^{ab}_e~r_a~s_b~v^e$ and $B^c_d~t_c~u^d$ and then multiplies those two together to get its final result. We use straight-up juxtaposition of the sub-tensors to denote this outer product.
The next thing we get is contraction, which requires an axiom: given an $[m, n]$-tensor there is some decomposition of it in terms of a sum of big outer products of $m$ vectors and $n$ covectors. Once this exists, you can just apply one of the covectors to one of the vectors, that creates a scalar, and a scalar times a tensor is just a tensor. As you can imagine, we denote this by repeating an index both-top-and-bottom. So $A^{mn}_m$ is a contraction which lives now in $\mathcal T^n$. It comes from some tensor $A^{an}_b$ by way of this axiom: that tensor was some sum of outer products: $$ A^{an}_b = \alpha^a~\beta^n~\gamma_b + \dots + \chi^a~\psi^n~\omega_b. $$ The resulting contraction is the vector $(\gamma_m~\alpha^m)\beta^n + \dots + (\omega_m~\chi^m)~\psi^n.$ Those terms in parentheses are applications of a covector to a vector to produce a scalar, so we're just adding together scalar multiples of vectors to create a new vector. So this is the coordinate-independent notion of trace which can be used to reduce any $[m, n]$-tensor to an $[m-1,n-1]$-tensor field, as long as neither of those two resulting numbers is negative.
Note: I took for granted above an important relabeling isomorphism, which identifies which $\alpha^m$ is the correct vector to represent $\alpha^a$ in a totally different space (one lives in $\mathcal T^a$, one lives in $\mathcal T^m$). We can write this relabeling isomorphism more explicitly as the tensor $\delta^a_m$ since it describes how to take any covector in $\mathcal T_a$ that would have operated directly on $\alpha^a$, and makes it operate instead on a vector in $\mathcal T^m.$ So, this is a coordinate-independent version of the "Kronecker delta."
Dot and cross products, $\nabla$
Finally dot products and cross products need to be implemented more directly. Well, that's easy: they are tensors in the space!
The dot product is a multilinear map from two vectors to a scalar, so it is a $[0, 2]$-tensor $g_{ab},$ called the metric tensor. It is symmetric in its two inputs, $g_{ab} = g_{ba}.$ It has an inverse $g^{ab}$ such that $g^{ab}~g_{bc} = \delta^a_c,$ and its inverse is naturally also symmetric. Just like the relabeling isomorphism is a canonical identification of different vector spaces, the metric is a canonical identification of vectors with covectors, every vector $\vec v$ corresponds to the covector $(\vec v\cdot)$ canonically. We usually indicate this by using the same symbol for the vector, e.g. $v_a = g_{ab}~v^b.$ We can do this with tensors too; usually we try to keep the indices in the same horizontal position while changing their vertical position from up to down, so that for example $$M^{ij}_{~~~k} = g^{ia}~g^{jb}~M_{ijk}.$$
An orientation tensor in $D$ dimensions is a totally antisymmetric $[0, D]$-tensor, so in 3 dimensions it is $\epsilon_{abc},$ implementing the cross product between two vectors. It is usually taken so that $$\epsilon^{ij\dots n}~\epsilon_{ij\dots n} = g^{ai}~g^{bj}~\dots g^{fn}~\epsilon_{ij\dots n}~\epsilon_{ab\dots f} = \pm D!,$$which usually makes the individual components $\pm 1$ in some suitable basis. The choice of $\pm$ in $\pm D!$ is usually related to the determinant of $g$ or some such, so when you get to Minkowski space I think -24 is more often used, but in 3D space +6 is used and in 4D Euclidean you'd probably use +24 instead.
For reasons that I can't explain due to answer length, $\nabla_i$ requires some more careful definition: the geometric structure it embodies is called a "connection" on the space, and it turns out that we have uniquely defined its operation on scalar fields, $v^i~\nabla_i~s$ being $V s$ directly: but its definition on vector fields is subject to an ambiguity up to a $[1, 2]$ tensor field (when you transport a vector along another vector, what vector results?). One particular one can be chosen, however, to have no "torsion" and to make the metric tensor have no derivative, this is called the "Levi-Civita connection," and this is usually what we use. Anyway, all I wanted to say here is that it clearly has a natural covector index, because of how we use it on scalars.
Introducing explicit coordinates again
All of the above indices are abstract, they just signify belonging to a set and allow some creative bookkeeping. But eventually you're going to want to introduce some sort of coordinates or so over the space and do some real calculations. When you want to do this, you basically construct within the space your local basis vectors $c_1^a, c_2^a, \dots c_D^a.$ It helps to choose two non-overlapping sets of symbols, e.g. Greek indices are abstract indices, Roman ones are stand-ins for numbers. So these coordinate vectors are $c_a^\alpha$.
We construct a set of "dual vectors" to these, really covectors, $c^a_\alpha$. The idea is that these need to be perpendicular to all of the ones that are not their "target" vector, and they need to be scaled so that their product with their target vector is 1: $$ c^m_{\alpha}~c_n^\alpha = \delta^m_n. $$ This is now your normal Kronecker $\delta$ symbol on the right, it's not some sort of clever isomorphism. Now any vector $v^\alpha$ can be operated on by these dual vectors to get some components $v^a = c^a_\alpha~v^\alpha$ which fully characterize it in terms of these fields. You can then construct the original from these components, $v^\alpha = \sum_a v^a~c_a^\alpha.$
By understanding the operations of $\nabla_\alpha$ on the component fields, one gets things like the Christoffel symbols.
Now in 3D flat Euclidean space, every scalar field is a smooth function of the $x, y, z$ fields of the space, and every vector field has canonical scalar field components $v_x, v_y, v_z$ defined over the whole space. If you write those as a column vector then covectors are row-vectors and our usual Euclidean metric identifies every vector with its covector as a transpose.
As a consequence $g_{mn} = c^\mu_m c^\nu_n g_{\mu\nu}$ for example is a Kronecker delta as well, and so one typically forgets entirely the distinction between lower and upper indices, the vector product $\vec u \cdot \vec v$ is mathematically equal to $\sum_i u_i~v_i,$ and we can even erase that $\sum_i$ if we use the Einstein summation convention. So you now know that in the "pure land" that these components come from, every lower index needs to be balanced by an upper one; but in this simpler world it is much easier, they just need to be paired to imply the summation over their components which entails a geometrically meaningful product. In this sense should you regard $\epsilon_{ijk}~\partial_j~B_k$: since $j$ and $k$ are repeated, this is implicitly summed over those components; the spare index $i$ indicates that what remains are the components of a $[1, 0]$-tensor.
So the "geometrically correct" version would be $\epsilon^{\lambda\mu}_{~~~~\nu}~\nabla_\mu~B^\nu,$ or so; but we know that we're in flat 3D space so we know that there is no harm from just dropping this pretense, as long as we get skittish around things like "I repeated this index three times, does that cause an Einstein summation?" (Answer: Yes, but WTF are you doing?)