Notation explanation for $\text{div} \ \omega = \nabla^i \omega_i $ and $\Delta = -\nabla^i \nabla_i $. 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.$
For a function, $\omega$, of n variables, $x_1$, $x_2$, ..., $x_n$, "grad $\omega$", or $\nabla \omega$, is defined as $\frac{\partial \omega}{\partial x_1}+ \frac{\partial \omega}{\partial x_2}+ \cdot\cdot\cdot+ \frac{\partial \omega}{\partial x_n}$.  That can be thought of as the "scalar product" of the "vector operator" $\nabla= \frac{\partial}{\partial x_1}+ \frac{\partial}{\partial x_2}+ \cdot\cdot\cdot+ \frac{\partial}{\partial x_n}$ with $\omega$ which is denoted, in Einstein summation notation, as "$\nabla_i \omega$".  

Similarly, divergence of a vector valued function (one-form) $\vec{\omega}= f_1(x_1,x_2,...,x_n)\vec{i_1}+ f_2(x_1,x_2,...,x_n)\vec{i_2}+ \cdot\cdot\cdot+ f_n(x_1,x_2,...,x_n)\vec{i_n}$, where "$i_m$" is the mth basis vector, is the scalar vector $\frac{\partial f_1}{\partial x_1}+ \frac{\partial f_2}{\partial x_2}+ \cdot\cdot\cdot+ \frac{\partial f_n}{\partial x_n}$ which can be thought of as the dot product of the vector operator $\nabla$ and the vector $\vec{\omega}$.  That is denoted, in Einstein summation notation, as $\nabla_i\omega^i$.