Sensitivity analysis is of significant importance for optimization, uncertainty quantification and mesh adaptation. While the conventional method of linearization and sensitivity analysis is well established for non-chaotic flows, it yields unbounded sensitivities for chaotic flows due to the presence of positive Lyapunov exponents [1]. Since the high-fidelity simulation of turbulent flows is inherently chaotic, it is of interest to obtain accurate sensitivities for such cases. In the past decade, shadowing-based approaches have been developed to compute accurate sensitivities for chaotic dynamical systems. Notable examples include the least-squares shadowing [2], non-intrusive least-squares shadowing [3], stabilized march [4], etc. The stabilized march approach computes the shadow trajectory by specifying suitable boundary conditions for the adjoint differential equation, and is provably convergent to the correct sensitivity with large integration times [4]. The current work uses stabilized march to compute the sensitivity for a three-dimensional test case involving turbulent flow. Comparisons against finite differences verify the accuracy of the approach. In addition, the convective term of the governing Navier-Stokes equation is discretized using the entropy-stable scheme, which discretely satisfies the entropy inequality [5]. An investigation is conducted on the impact of entropy-stable schemes, which are known to lack the local energy stability property [6], on the sensitivity of chaotic flows.