**Strategy:** We proceed by a density argument. Since the boundary is $C^1$, functions in $C^1(\bar{U})$ are dense in $W^{1,p}(U)$. We will define the operator $E$ for smooth functions, prove the explicit norm estimates, and then extend $E$ to all of $W^{1,p}(U)$ by the Bounded Linear Transformation theorem.

**Step 1: Local Boundary Flattening** Since $\partial U$ is compact and of class $C^1$, there exists a finite open cover $\{W_i\}_{i=1}^N$ of $\partial U$ and corresponding $C^1$-diffeomorphisms (charts). Let $W_0 \subset \subset U$ be an open set such that $U \subset W_0 \cup \bigcup_{i=1}^N W_i$. For each $i \in \{1, \dots, N\}$, let $\Phi_i: W_i \to Q$ be a $C^1$ diffeomorphism mapping $W_i$ onto the unit cube $Q = (-1, 1)^n$, such that: \begin{align*} \Phi_i(W_i \cap U) &= Q_+ := \{ y \in Q : y_n > 0 \} \\ \Phi_i(W_i \cap \partial U) &= Q_0 := \{ y \in Q : y_n = 0 \}. \end{align*} Let $\{\zeta_i\}_{i=0}^N$ be a smooth partition of unity subordinate to the cover $\{W_i\}$.

**Step 2: Higher-Order Reflection** We define the extension in the flattened coordinates. Let $v \in C^1(\bar{Q}_+)$. We extend $v$ to $Q_- := \{y \in Q : y_n < 0\}$ via a higher-order reflection $E_0 v$. Define the reflection map: \begin{align*}

E_0 v(y', y_n) := \begin{cases} v(y', y_n) & \text{if } y_n \ge 0 \\ -3v(y', -y_n) + 4v(y', -y_n/2) & \text{if } y_n < 0. \end{cases}

\end{align*}

**Step 3: Explicit Norm Estimates for Reflection** We verify the $W^{1,p}$ boundedness of $E_0$. First, the $L^p$ norm on the lower half $Q_-$. Using the change of variables $z_n = -y_n$ (Jacobian determinant is 1) and $w_n = -y_n/2$ (Jacobian determinant is 2), we compute: \begin{align*} \int_{Q_-} |E_0 v|^p \, \mathrm{d}\mathcal{L}^n(y) &= \int_{Q_-} |-3v(y', -y_n) + 4v(y', -y_n/2)|^p \, \mathrm{d}\mathcal{L}^n(y) \\ &\le 2^{p-1} \left( 3^p \int_{Q_-} |v(y', -y_n)|^p \, \mathrm{d}y + 4^p \int_{Q_-} |v(y', -y_n/2)|^p \, \mathrm{d}y \right) \\ &= 2^{p-1} \left( 3^p \int_{Q_+} |v(y', z_n)|^p \, \mathrm{d}z + 4^p \int_{Q_+ \cap \{w_n < 1/2\}} |v(y', w_n)|^p \cdot 2 \, \mathrm{d}w \right) \\ &\le C_0 \|v\|_{L^p(Q_+)}^p. \end{align*} Next, we verify the derivatives. For tangential derivatives $j < n$: \begin{align*} \partial_{y_j} (E_0 v) = -3 (\partial_{y_j} v)(y', -y_n) + 4 (\partial_{y_j} v)(y', -y_n/2). \end{align*} This has the exact same structure as the function values, so the $L^p$ norm is bounded by $C \|\partial_{y_j} v\|_{L^p(Q_+)}$. For the normal derivative $j=n$, using the Chain Rule: \begin{align*} \partial_{y_n} (E_0 v) &= -3 (\partial_{y_n} v)(y', -y_n) \cdot (-1) + 4 (\partial_{y_n} v)(y', -y_n/2) \cdot (-1/2) \\ &= 3 (\partial_{y_n} v)(y', -y_n) - 2 (\partial_{y_n} v)(y', -y_n/2). \end{align*} Again, performing the same change of variables ($z_n = -y_n, w_n = -y_n/2$): \begin{align*} \|\partial_{y_n} (E_0 v)\|_{L^p(Q_-)}^p \le C \left( \|\partial_{y_n} v\|_{L^p(Q_+)}^p + \|\partial_{y_n} v\|_{L^p(Q_+)}^p \right). \end{align*} Summing these estimates: \begin{align*} \|E_0 v\|_{W^{1,p}(Q)} \le C \|v\|_{W^{1,p}(Q_+)}. \end{align*}

**Step 4: Global Synthesis and Stability** Let $u \in C^1(\bar{U})$. We decompose $u$ using the partition of unity: $u = \sum_{i=0}^N \zeta_i u$. 1. **Interior Term ($i=0$):** The function $\zeta_0 u$ has support in $W_0 \subset U$. We extend it by 0 to $\mathbb{R}^n$. \begin{align*} \|\zeta_0 u\|_{W^{1,p}(\mathbb{R}^n)} = \|\zeta_0 u\|_{W^{1,p}(U)} \le C_{\zeta} \|u\|_{W^{1,p}(U)}. \end{align*} 2. **Boundary Terms ($i \ge 1$):** For each patch $W_i$, define $u_i := \zeta_i u$. We define the extension $Eu_i$ by mapping to $Q$, reflecting, and mapping back: \begin{align*} Eu_i := (E_0(u_i \circ \Phi_i^{-1})) \circ \Phi_i. \end{align*} Since $\Phi_i$ and $\Phi_i^{-1}$ are $C^1$ on compact sets, their Jacobians and derivatives are uniformly bounded. Let $K$ be a bound on $\|\Phi_i\|_{C^1}$ and $\|\Phi_i^{-1}\|_{C^1}$. Using the Chain Rule and Change of Variables formula for Sobolev functions: \begin{align*} \|Eu_i\|_{W^{1,p}(W_i)} &\le K_1 \|E_0(u_i \circ \Phi_i^{-1})\|_{W^{1,p}(Q)} \\ &\le K_1 C_{\text{refl}} \|u_i \circ \Phi_i^{-1}\|_{W^{1,p}(Q_+)} \\ &\le K_1 C_{\text{refl}} K_2 \|u_i\|_{W^{1,p}(W_i \cap U)}. \end{align*} Thus, $\|Eu_i\|_{W^{1,p}(\mathbb{R}^n)} \le C_i \|u\|_{W^{1,p}(U)}$.

**Step 5: Definition of $E$** Define the total extension for smooth $u$: \begin{align*} Eu := (\zeta_0 u)_{\text{zero ext}} + \sum_{i=1}^N Eu_i. \end{align*} By the triangle inequality and the estimates in Step 4: \begin{align*} \|Eu\|_{W^{1,p}(\mathbb{R}^n)} \le \sum_{i=0}^N C_i \|u\|_{W^{1,p}(U)} = C_{\text{total}} \|u\|_{W^{1,p}(U)}. \end{align*} Since $C^1(\bar{U})$ is dense in $W^{1,p}(U)$, the operator $E$ extends uniquely to a bounded operator on all of $W^{1,p}(U)$.