Divergence of a function product with functions in $L^infty (Omega)$ and $H^1(Omega)$

Let me reach a bit further back.

Your notation is a bit sloppy, which may mislead and thus confuse you. The correct expression for the differential operator is $nabla cdot (k nabla u) = operatorname{div} (k nabla u)$, i.e. the divergence of $k nabla u$, while $nabla (k nabla u)$ would denote the gradient of the vector field $k nabla u$, i.e. the derivative in all directions (the Jacobian). This distinction is important because requiring $operatorname{div}(k nabla u)$ to be well-defined is a much weaker condition than requiring $nabla(k nabla u)$ to be well-defined. In the following I assume you wanted to ask what conditions are necessary for $operatorname{div}(k nabla u)$ to be well-defined.

To be nit-picky (or maybe not), you must specify in which sense that expression should be well-defined. I assume you meant that $operatorname{div}(k nabla u) in L^2(Omega)$ since then, under the additional assumption that $f in L^2(Omega)$, one can say that the differential equation $operatorname{div}(k nabla u(x)) = f(x)$ holds for $x in Omega$ almost everywhere.

If you don't assume from the get-go that $k nabla u in (H^1(Omega))^d$ then the question arises what $operatorname{div}(k nabla u) in L^2(Omega)$ should even mean. You stated in your comments that you never found a use for distributional derivatives, so this is a good time as any to make use of those concepts (and maybe clear up some more confusion). For any vector field $v in (L^2(Omega))^d$, with the special case $v = k nabla u$ in mind, the distributional divergence of $v$ is the distribution $operatorname{div} v in C^infty_0(Omega)'$ defined by $$ [operatorname{div} v](varphi) := -int_{Omega} v cdot nabla varphi , mathrm dx quad forall varphi in C^infty_0(Omega). $$

This mirrors the definition of the distributional derivatives which you should already know. Indeed, you can define the distributional version of any linear differential operator in a similar way. It's important to realize that while this may look like "just" partial integration / Green's identities, the basic idea of distributional derivatives is that the right-hand side is well-defined even if only $v in (L^2(Omega))^d$ and nothing more. Only if $v$ is more regular, partial integration allows us to a-posteriori identify the distributional derivative with more a classical function. This identification of distributional divergence is done (as for distributional derivatives) by saying that $operatorname{div} v in C^infty_0(Omega)'$ is represented by a function $h in L^2(Omega)$ if there holds $$ int_{Omega} h varphi ,mathrm dx = v(varphi) left( = -int_{Omega} v cdot nabla varphi , mathrm dx right)$$ for all $varphi in C^infty_0(Omega)$. One then, casually one may say, writes $operatorname{div} v = h$. This is the sense in which $operatorname{div} v in L^2(Omega)$, e.g. $operatorname{div}(k nabla u) in L^2(Omega)$, must be understood. Of course one could just say that a vector field $v$ has "a divergence in $L^2(Omega)$" if there is a $h in L^2(Omega)$ satisfying this equation without ever mentioning "distribution" or "distributional derivative", but this hides the generality of the ideas behind the construction.

With this preparation, I claim that $k in L^infty(Omega)$ is sufficient that any solution $u in H^1(Omega)$ of your weak equation $$ int_Omega k nabla u cdot nabla varphi ,mathrm dx = int_Omega f varphi ,mathrm dx $$ for all $varphi in H^1_0(Omega)$ with $f in L^2(Omega)$ has regularity $operatorname{div}(k nabla u) in L^2(Omega)$. Comparing the integral equation for the weak solution with the condition I wrote above for $h in L^2(Omega)$ to represent the distributional derivative of $operatorname{div}(k nabla u)$, you can see that $h = - f$ and if $f in L^2(Omega)$ thus $operatorname{div}(k nabla u) in L^2(Omega)$. The regularity $k in L^infty(Omega)$ is here only used such that $v = k nabla u in L^2(Omega)$ (and of course for less regular $k$ the integral equation for the weak solutions is no longer well-defined for $varphi in H^1_0(Omega))$). This argument of course only works because $u$ is not any function in $H^1(Omega)$ but one that satisfies the integral equation, which gives you the regularity $operatorname{div}(k nabla u) in L^2(Omega)$ "for free" if your right-hand side $f$ has that regularity.

As a last note, the space of vector fields $v in (L^2(Omega))^d$ with $operatorname{div} v in L^2(Omega)$ is typically denoted by $H(operatorname{div}; Omega)$, thus $k nabla u in H(operatorname{div}; Omega)$ is the minimal regularity for $operatorname{div}(k nabla u)$ to be well-defined.

I believe you would need $k nabla u in H^1$ to have these defined directly. Assuming that $k in L^infty$ and $nabla u in H^1$ in general isn't enough to make the product in $H^1$ (take $k$ to be very not-differentiable, but still bounded, and take $nabla u$ to be nonzero where $k$ is badly behaved). However, you could assume the weaker condition that $k nabla u in L^2$, and then define $nabla (k nabla u)$ weakly by the formula $$ int_Omega phi nabla(k nabla u), dx = -int_Omega nabla phi cdot k nabla u, dx $$ for all $phi$ smooth and compactly supported in $Omega$ (or all $phi in H^1_0(Omega))$. The condition $knabla u in L^2$, for instance, would follow if $u in H^1$ and $k in L^infty$.

When doing weak formulation of PDE, it is usually this second approach that is used. The unknown $u$ is generally required to be $H^1(Omega)$ (or $H^1_0(Omega)$ if we have Dirichlet boundary conditions), and the coefficients $k$ may by $L^infty$ or whatever we suppose for the problem.