Embedded boundary (EB) methods provide a sharp representation of solid boundaries by cutting out the solid body from a Cartesian fluid grid. They highly simplify the mesh generation for flow simulations over practical geometries and flow conditions. Turbulent flows in supersonic regimes involve broadband flow scales (that are highly sensitive to numerical dissipation) as well as strong shocks (that may generate spurious oscillations with non-dissipative schemes). Robust computations of such flows require upwind (or biased) schemes in flow regions with discontinuities and non-dissipative centered schemes in the smooth flow regions. Weighted essentially non-oscillatory (WENO) schemes are commonly used to adapt the stencils based on the smoothness of flow variables. In bounded domains, the WENO schemes are often supplemented by low-order boundary closures or ghost point treatments that can compromise the accuracy and the conservation properties of the overall scheme. At embedded boundaries, an ad hoc boundary treatment can result in numerical instabilities detrimental to long-time flow calculations. This study uses energy stability concepts to develop provably stable EB closures for a sixth-order WENO scheme. In the proposed finite-difference discretization, the solutions are stored and advanced at the grid points (or the cell centers), while the fluxes are computed at cell interfaces (or the flux points) to ensure conservation by construction. The smoothness indicator is derived to obtain an energy estimate for the overall scheme. The small-cell problem that arises when the cut cell adjoining the EB is significantly smaller than a regular cell is addressed by choosing the flux point spacings that do not vanish even when the grid points at the EB coincide. This procedure avoids additional steps, e.g. cell mixing/merging, flux redistribution, etc., that are difficult to automate. The derived scheme is dimensionally-split and, hence, does not require geometry/solution reconstructions, making it easy to incorporate in an existing solver. Various linear and non-linear, inviscid and viscous numerical tests are performed to demonstrate the stability and the accuracy of the proposed method.