Given a connection $\nabla$ on the tangent bundle, the divergence of a vector field $X\in C^{\infty}(TM)$ is defined as $div(X) = C_1^1(\nabla X)$ where $C$ is the contraction.
If $\{e_i\}$ is a local reference frame and $\{\eta^i\}$ is the induced coframe, then locally $$(1)\quad \quad div(X) = \nabla X(e_i,\eta^i)=\eta^i(\nabla_i X) = e_i(X^i) + X^r\Gamma_{i r}^i $$ where $\Gamma_{i r}^k$ are the Christoffel's symbols.
Given a pseudo-Riemannian metric, if instead of a vector field we have a 1-form $\omega\in C^{\infty}(T^*M)$, we can define $$(2) \quad \quad div(\omega) := div(\omega^\#).$$
Sometimes I have encountered this notation for the divergence of a 1-form : $$ (3)\quad \quad \text{div} \ \omega = \nabla^i \omega_i$$ What does the symbol $\nabla^i$ mean? Is $\omega_i$ defined by $\omega = \omega_i dx^i$? Or it is $\text{div} \ \omega = (\nabla^i \omega)_i$ the i-component of the covariant tensor $\nabla^i \omega$?
I also need to explain why , in the same paper, the author states that $\Delta = -\nabla^i\nabla_i$ as if $\Delta f = -\nabla^i\nabla_if$ since $\nabla_i f $ is a function it seems that $\nabla^i$ act on functions.
Locally $(2)$ is written as $div(\omega) = e_i( g^{ij}\omega_j)+ g^{r j}\omega_j\Gamma_{i r}^i $ I can't see clearly how is defined $\nabla^i$.
This is index notation as traditionally used in the Ricci calculus, these days typically interpreted by mathematicians as abstract index notation.
When we write $\mathrm{div}\,X = \nabla_i X^i$, what is meant is that we first take the covariant derivative of $X$ and then perform the contraction indicated by the repeated indices; i.e. $\nabla_i X^i = (\nabla X)^i_i$, where the RHS can be interpreted either as abstract index notation denoting the contraction $C(\nabla X)$ or as a coordinate formula using Einstein summation notation. Importantly, this is not the same thing as first taking the coordinate components $X^i$ (which are scalars) and then differentiating, which would yield instead $\partial_i X^i$, a coordinate-dependent quantity. Thus we usually don't think of the expression $\nabla_i X^i$ as being an operator $\nabla_i$ acting on the component functions $X^i$ - instead, it's the operator $\nabla$ acting on the vector field $X$, with the indices just telling us how to contract. (Hopefully this answers your last question regarding $\nabla^i \nabla_i f.$)
In the case of a one-form $\omega$, the covariant derivative is a $(0,2)$-tensor $\nabla\omega = (\nabla_j \omega_i) \eta^j \otimes \eta^i$. (Note again that in this notation, $\nabla_j \omega_i$ is not the same as $\partial_j \omega_i$ - we are taking the covariant derivative before plugging in the reference frame.) To take the trace of this to obtain a divergence, we first raise an index with the metric as you described, yielding $\mathrm{div}(\omega^\sharp) = \nabla_i \omega^i := \nabla_i (\omega^\sharp)^i$ where of course $\omega^i = (\omega^\sharp)^i = g^{ij} \omega_j.$ Since the covariant derivative is $g$-compatible and the metric is symmetric, we can think of this as $\mathrm{tr}_g \nabla \omega = g^{ij} (\nabla \omega)_{ij};$ so we could just as fairly raise the first index instead of the second, yielding $\nabla^i \omega_i.$