A direct sensitivity analysis technique is extended to calculate higher-order sensitivity coefficients in three-dimensional air quality models. The time evolution of sensitivity coefficients of different order is followed alongside that of the concentrations. Calculation of higher-order sensitivity coefficients requires few modifications to the original (first-order) sensitivity modules and is carried out efficiently and with minimal computational overhead. The modeling results (first-, second-, and third-order sensitivity coefficients) for an ozone episode in central California are shown and discussed. Second-order sensitivity coefficients of ozone concentration with respect to domain-wide NO emissions show reasonable agreement with brute-force results and exhibit less noisy behavior. By using second-order sensitivity coefficients the nonlinear responses are better captured and described. For a Taylor series projection from the base case, including the second-order term improves the accuracy. In general, higher-order sensitivity analysis shows a noticeable improvement in terms of accuracy over the conventional first-order analysis. Of particular interest, second-order sensitivity analysis is better equipped to address the nonlinear behavior around the peak ozone in NO(x)-rich plumes.