In the numerical gray zone of the convective boundary layer (CBL), the horizontal resolution is comparable to the size of organized convective circulation. As turbulence becomes partially resolved, grid-scale variations of the subgrid-scale (SGS) turbulent fluxes become significant compared to the mean. Previously, such variations have often been ignored in scale-adaptive planetary boundary layer (PBL) schemes developed for the gray zone. This study investigates these variations with respect to height and resolution based on large-eddy simulations. It is found that SGS fluxes exhibit maximum variability at the center of the gray zone, where the resolved and the SGS mean fluxes are approximately equal. A simple analytical model is used to associate such characteristic variations to the nonlinear interactions of the dominant energy-containing mode of CBL turbulence. Examination of the horizontal distribution of the SGS fluxes reveals their preferential location over the updraft edges surrounding the core. A priori analysis further suggests the ability of a scale-similarity closure to reproduce the unique spatial patterns of the SGS fluxes at gray zone resolutions. Four scale-adaptive PBL schemes are evaluated focusing on their representations of the modeled SGS flux variability. Their shared short-comings as a result of their gradient diffusion-based formulation are exposed. This study suggests that a mixed model consisting of a scale-adaptive PBL scheme to represent the mean, and a scale-similarity component to account for grid-scale variability to be advantageous for the gray zone.