For a multilinear mapping, it suffices to consider its Frechet derivative. Let $W$ be an $n$-D vector space, and each $V_i$ be an $m_i$-D vector space with $i=1,2,...,N$. Let $f:V_1\times V_2\times\cdots\times V_N\to W$ be multilinear. Then $\forall\left(v_1,v_2,...,v_N\right)\in V_1\times V_2\times\cdots\times V_N$, the Frechet derivative of $f$ at this location, denoted by $({\rm d}f)(v_1,v_2,...,v_N)$, is also a multilinear mapping, i.e.,
$$
({\rm d}f)(v_1,v_2,...,v_N):V_1\times V_2\times\cdots\times V_N\to W.
$$
According to Frechet, it follows that
\begin{align}
&({\rm d}f)(v_1,v_2,...,v_N)(h_1,h_2,...,h_N)\\
&=f(h_1,a_2,...,a_N)\\
&+f(a_1,h_2,...,a_N)\\
&+\cdots\\
&+f(a_1,a_2,...,h_N).
\end{align}
Recall that, if $g$ is linear, its entry-wise form reads
$$
g_i(v)=\sum_ja_{ij}v_j,
$$
and if $g$ is bilinear, its entry-wise form reads
$$
g_i(v_1,v_2)=\sum_{j_1,j_2}a_{ij_1j_2}v_{1j_1}v_{2j_2}.
$$
Inductively and formally, the above multilinear $f$ observes the following entry-wise form
$$
f_i(v_1,v_2,...,v_N)=\sum_{j_1=1}^{m_1}\sum_{j_2=1}^{m_2}\cdots\sum_{j_N=1}^{m_N}a_{ij_1j_2...j_N}v_{1j_1}v_{2j_2}...v_{Nj_N}
$$
for $i=1,2,...,m$, where each $v_{kj_k}$ denotes the $j_k$-th entry of $v_k\in V_k$, while $a_{ij_1j_2...j_N}$'s are the coefficients of $f$.
Thanks to this entry-wise form, we may then write down the entry-wise form of ${\rm d}f$ as well, which reads
\begin{align}
&({\rm d}f)_i(v_1,v_2,...,v_N)(h_1,h_2,...,h_N)\\
&=\sum_{j_1=1}^{m_1}\sum_{j_2=1}^{m_2}\cdots\sum_{j_N=1}^{m_N}a_{ij_1j_2...j_N}h_{1j_1}v_{2j_2}...v_{Nj_N}\\
&+\sum_{j_1=1}^{m_1}\sum_{j_2=1}^{m_2}\cdots\sum_{j_N=1}^{m_N}a_{ij_1j_2...j_N}v_{1j_1}h_{2j_2}...v_{Nj_N}\\
&+\cdots\\
&+\sum_{j_1=1}^{m_1}\sum_{j_2=1}^{m_2}\cdots\sum_{j_N=1}^{m_N}a_{ij_1j_2...j_N}v_{1j_1}v_{2j_2}...h_{Nj_N}.
\end{align}
In other words, as $a_{ij_1j_2...j_N}$'s are known, the entry-wise form of ${\rm d}f$ could be expressed straightforwardly as above.
Finally, the "$+$" in OP's original post, i.e., $(h_1+h_2+\cdots+h_N)$, is a convention in some context, which is exactly $(h_1,h_2,...,h_N)$ here. When there is free of ambiguity, both expressions can be used as per ones preference.
Best Answer
Defining a "del" vector $\nabla=\mathbf{e}_i\frac{\partial}{\partial x^i}$ (summation on repeated indices assumed) is sort-of a mistake. This isn't a vector (field). It is just a nice notational device, because in this case, in cartesian coordinates, we have $$ \text{grad}(f)=\mathbf{e}_i\frac{\partial f}{\partial x^ i},$$ which looks like $\nabla f$, where $\nabla$ is the "vector" defined above and we have $$ \text{div}(\mathbf{A})=\frac{\partial A^ i}{\partial x^i}, $$ which sort of looks like $\nabla\cdot\mathbf{A}$, etc. But this is just a useful way to remember these formulas.
Now let's note several things. We can define a matrix whose $i,j$-th element is $\partial A^ i/\partial x^ j$. Let's call this matrix $\nabla\otimes\mathbf{A}$. We clearly have $$ \text{div}(\mathbf{A})=\text{Tr}(\nabla\otimes\mathbf{A})=\sum_i\frac{\partial A^ i}{\partial x^ i}. $$ (I'm gonna suspend the automatic summation from now on.)
If we have an orthogonal curvilinear coordinate system $(u^1,u^2,u^3)$, we can define the following quantities:
Coordinate basis vectors $$ \mathbf{g}_i=\frac{\partial \mathbf{r}}{\partial u^ i}, $$
"Reciprocal" coordinate basis vectors $$ \mathbf{g}^i=\nabla u^ i, $$
An orthonormal frame $$ \hat{\mathbf{e}}_i=\frac{1}{h_i}\mathbf{g}_i, $$ where $h_i=\sqrt{\mathbf{g}_i\cdot\mathbf{g}_i}$.
These quantities have the property that $$ \mathbf{g}^ i\cdot\mathbf{g}_j=\delta^ i_j, $$ since ($x^ i$ are cartesian coordinates and $\mathbf{e}_i$ are cartesian basis vectors) we have $$ \mathbf{g}^i\cdot\mathbf{g}_j=\nabla u^i\cdot\frac{\partial\mathbf{r}}{\partial u^ j}=\sum_k\frac{\partial u^ i}{\partial x^k}\mathbf{e}_k\cdot\sum_l\frac{\partial x^l}{\partial u^ j}\mathbf{e}_l=\sum_{kl}\frac{\partial u^ i}{\partial x^k}\frac{\partial x^l}{\partial u_j}\delta_{kl}=\sum_k \frac{\partial u^ i}{\partial x^k}\frac{\partial x^k}{\partial u_j}=\delta^i_j.$$
We further introduce a matrix $g_{ij}$ given as $g_{ij}=\mathbf{g}_i\cdot\mathbf{g}_j$. Because of orthogonality, we have $g_{ij}=h_i^2\delta_{ij}$. We also introduce $g^{ij}=\mathbf{g}^i\cdot\mathbf{g}^j$, feel free to check that $g^{ij}=\frac{1}{h_i^2}\delta_{ij}$, thus $g_{ij}$ and $g^{ij}$ are inverse matrices.
If $\mathbf{A}$ is a vector field, we associate three kinds of components with $\mathbf{A}$. We can have
$\mathbf{A}=\sum_i A^ i\mathbf{g}_i$, and the $A^ i$ are called contravariant components,
$\mathbf{A}=\sum_i A_i\mathbf{g}^i$, and the $A_i$ are called covariant components and
$\mathbf{A}=\sum_i \hat{A}_i\hat{\mathbf{e}}_i$ and the $\hat{A}_i$ are called physical components (here whether you put the index upwards or downwards doesn't matter!).
It is clear that we have $A^ i=\mathbf{g}^i\cdot\mathbf{A}$ and $A_i=\mathbf{g}_i\cdot\mathbf{A}$, and what is less clear is that differentiation with respect to the coordinates $\partial/\partial u^ i$ produce covariant components, eg. we have for any scalar function $f$, $\nabla f=\sum_i\mathbf{g}^i\frac{\partial f}{\partial u^ i}$ (check the definition of $\mathbf{g}^i$!).
We define Christoffel symbols of the second kind as $$ \Gamma^ k_{ij}=\frac{1}{2}\sum_l g^{kl}\left(\frac{\partial g_{jl}}{\partial u^ i}+\frac{\partial g_{il}}{\partial u^ j}-\frac{\partial g_{ij}}{\partial u^ l}\right). $$ Check that $\Gamma^k_{ij}=\mathbf{g}^k\cdot\frac{\partial \mathbf{g}_i}{\partial u^ j}$ (I really don't want to derive this identity right now).
In curvilinear coordinates, the basis vectors also depend on positions, so every time you differentiate a vector field, you need to make sure to take the variation of the basis vectors also into account, so we calculate the divergence as $$ \text{div}(\mathbf{A})=\sum_i\mathbf{g}^i\cdot\frac{\partial\mathbf{A}}{\partial u^i}=\sum_{ij}\mathbf{g}^i\left(\frac{\partial A^j}{\partial u^ i}\mathbf{g}_j+A^j\frac{\partial \mathbf{g}_j}{\partial u^ i}\right)=\sum_{ij}\left(\frac{\partial A^ i}{\partial u^i}+\Gamma^i_{ji}A^j\right). $$
We now calculate $$\Gamma^i_{ji}=\sum\frac{1}{2}g^{il}\left(\frac{\partial g_{jl}}{\partial u^i}+\frac{\partial g_{il}}{\partial u^j}-\frac{\partial g_{ij}}{\partial u^l}\right)=\sum\frac{1}{2}g^{il}\frac{\partial g_{il}}{\partial u^j},$$ which, in matrix notation ($(g_{ij})=g$) can be written as $$ \Gamma^i_{ji}=\frac{1}{2}\text{Tr}\left(g^{-1}\frac{\partial g}{\partial u^j}\right), $$ which, using Jacobi's formula, can be written as $$ \sum_i\Gamma^ i_{ji}=\frac{1}{2}\frac{1}{\det g}\frac{\partial}{\partial u^ j}\det g=\frac{1}{2}\frac{\partial}{\partial u^ j}\ln\det g=\frac{\partial}{\partial u^j}\ln\sqrt{\det g}=\frac{1}{\sqrt{\det g}}\frac{\partial}{\partial u^j}\sqrt{\det g}. $$
Plugging this back into the divergence formula gives $$\text{div}(\mathbf{A})=\sum_{ij}\left(\frac{\partial A^ i}{\partial u^i}+\Gamma^i_{ji}A^j\right)=\sum_{ij}\left(\frac{\partial A^i}{\partial u^i}+\frac{1}{\sqrt{\det g}}\frac{\partial}{\partial u^j}\sqrt{\det g}A^j\right)= \\ =\sum_{i}\frac{1}{\sqrt{\det g}}\frac{\partial}{\partial u^i}(\sqrt{\det g}A^i). $$
But since $g=\text{diag}(h_1^2,h_2^2,h_3^2)$
, we have $\det g=(h_1h_2h_3)^2$, so $$ \text{div}(\mathbf{A})=\frac{1}{h_1h_2h_3}\left\{\frac{\partial}{\partial u^1}(h_1h_2h_3A^1)+\frac{\partial}{\partial u^2}(h_1h_2h_3A^2)+\frac{\partial}{\partial u^3}(h_1h_2h_3A^3)\right\}= \\ =\frac{1}{h_1h_2h_3}\left\{\frac{\partial}{\partial u^1}(h_2h_3\hat{A}_1)+\frac{\partial}{\partial u^2}(h_1h_3\hat{A}_2)+\frac{\partial}{\partial u^3}(h_1h_2\hat{A}_3)\right\}, $$
which is the formula for divergence in terms of the physical components.