Part 2:
As pointed out in the answers, there is a paper by Doyoon Kim which completely answers the second part of my question. An earlier paper by Jürgen Marschall also answers this question, with possible exception of some exponents (interestingly, this paper appeared before Costabel's paper). I say "possible", because the both papers treat the traces for $L^p$-based spaces and Marschall's paper exclude the case $s-\frac1p$ is an integer for the boundedness of the right inverse.
But his methods maybe applicable for $p=2$ (I have not checked) as we know this case is somewhat special. In any case Kim proves the theorem even when $s-\frac1p$ is an integer, for general $p$.
Part 1:
The first part of the question, that is the claim about the endpoint case $s=k+\frac32$ is not true, in the sense that the norm of the trace operator cannot be bounded in terms of the Lipschitz constant alone. Some years after the paper containing the claim appeared, Jerison and Kenig corrected it by publishing a counterexample for domains with $C^1$ boundary, which means there is no hope for the Lipschitz case ($k=0$). Note that the main results in the earlier paper by Jerison and Kenig depend weakly on this claim, and are still true. Their $C^1$ counterexample was based on a counterexample in the Lipschitz case due to Guy David, and I want to record its main idea here.
Consider in $\mathbb{R}^2$ a half disk whose flat part of the boundary is replaced by a sawtooth (or zigzag) with period $2\varepsilon$.
Let us call this domain $\Omega_\varepsilon$.
Note that the Lipschitz constant of the domain does not depend on $\varepsilon$.
We pick a function $f\in H^{3/2}(\mathbb{R})$ such that $f'(y)$ tends to $\infty$ as $y\to0$. This is possible because $H^{1/2}(\mathbb{R})$ has unbounded functions. Then in the configuration as shown in the picture, we define $g(x,y)=\phi(x,y)f(y)$, where $\phi$ is a suitable cut-off function (say, equal to 1 on a considerable portion of the zigzag boundary).
Obviously, we have $g\in H^{3/2}(\Omega_\varepsilon)$,
but the tangential derivative of $g$ along the zigzag boundary is $\pm\sqrt2f'$ on a large portion of the boundary. This means that the $L^2$-norm of the tangential derivative is at least a constant multiple of $|f'(\varepsilon)|$, which blows up as $\varepsilon\to0$.
