This article introduces a novel approach for two-dimensional (2D) direction-of-arrival (DOA) estimation using stacked intelligent metasurfaces (SIM). SIMs are advanced structures composed of multiple metasurface layers, each containing numerous meta-atoms that can manipulate electromagnetic (EM) waves. The key innovation is the use of a SIM to perform a 2D discrete Fourier transform (DFT) as EM waves propagate through it, allowing the receiver array to directly observe the angular spectrum of the incoming signal. This eliminates the need for complex digital signal processing and reduces hardware complexity by avoiding power-hungry RF chains. The authors formulate an optimization problem to minimize the fitting error between the SIM's EM response and a 2D DFT matrix, and they develop a gradient descent algorithm to iteratively update the phase shifts of the meta-atoms. Additionally, the phase shift pattern in the zeroth layer is configured to generate orthogonal spatial frequency bins, improving DOA estimation accuracy. Theoretical analysis provides an upper bound for the mean square error (MSE) of the DOA estimator, and numerical simulations confirm that a well-trained SIM achieves an MSE of $10^{-4}$ for DOA estimation. The paper also discusses the system model, SIM optimization for 2D DFT, and the DOA estimation protocol. The proposed method offers a significant advantage over conventional DOA estimation techniques by leveraging the inherent parallel processing capabilities of SIMs, enabling fast and energy-efficient signal processing. The results demonstrate the potential of SIMs for future high-speed, low-power applications in 6G and beyond.This article introduces a novel approach for two-dimensional (2D) direction-of-arrival (DOA) estimation using stacked intelligent metasurfaces (SIM). SIMs are advanced structures composed of multiple metasurface layers, each containing numerous meta-atoms that can manipulate electromagnetic (EM) waves. The key innovation is the use of a SIM to perform a 2D discrete Fourier transform (DFT) as EM waves propagate through it, allowing the receiver array to directly observe the angular spectrum of the incoming signal. This eliminates the need for complex digital signal processing and reduces hardware complexity by avoiding power-hungry RF chains. The authors formulate an optimization problem to minimize the fitting error between the SIM's EM response and a 2D DFT matrix, and they develop a gradient descent algorithm to iteratively update the phase shifts of the meta-atoms. Additionally, the phase shift pattern in the zeroth layer is configured to generate orthogonal spatial frequency bins, improving DOA estimation accuracy. Theoretical analysis provides an upper bound for the mean square error (MSE) of the DOA estimator, and numerical simulations confirm that a well-trained SIM achieves an MSE of $10^{-4}$ for DOA estimation. The paper also discusses the system model, SIM optimization for 2D DFT, and the DOA estimation protocol. The proposed method offers a significant advantage over conventional DOA estimation techniques by leveraging the inherent parallel processing capabilities of SIMs, enabling fast and energy-efficient signal processing. The results demonstrate the potential of SIMs for future high-speed, low-power applications in 6G and beyond.