In this contribution we study model order reduction of nonlinear transport networks for district heating systems. In these, heating energy injected at a centralized power plant is transported to consumers. They are modeled by a hyperbolic differential algebraic system with large state space dimension. The network structure introduces sparse system dynamics, which transform to a dense reduced system leading to unacceptable computational costs [1]. To exploit the benefits of sparsity, sub-parts of the network are reduced separately in a structure preserving way using Galerkin projection [2]. After introducing the dynamical model, we discuss a decomposition strategy which aims at minimizing the number of observables for each subnetwork. We demonstrate its numerical benefits at an existing large scale heating network.The dynamics of heating networks are described by one dimensional incompressible Euler equations. The latter incorporate the advection of the energy density ϕ by a purely time dependent volume flow q(t) in each pipeline, as well as the conservation of momentum. In this, the pressure difference along a pipeline is dominated by frictional forces, allowing to neglect the contribution of acceleration in a quasi-stationary limit. After spatial discretization of the advection equation using the upwind scheme, a state space decomposition is performed allowing to formulate the differential algebraic equationThe algebraic constraints (3) result from the coupling of pipelines at nodes of the network and gather three contributions: the conservation of volume, Kirchhoff's second law claiming that pressure differences along network loops sum up to zero, and the volumeflow balance at consumers. It forces the product of volume flow q at consumer stations and observed energy density y to equal the power consumption u H . Hence, (3) collects quadratic constraints in the state variables ϕ, q. The system description is decomposed to a main network ϕ 0 to which the thermal control u T is applied, and c ∈ N subnetworks forming ϕ which receive observablesỹ as artificial inputs. Observables y are measured by both main-and subnetworks, whileỹ is the system state of the main network at attachment points to subnetworks.Ã,B,C represent continued block diagonal matrices. Their blocks A s , B s , C s , s ∈ {1, .., c} describe the dynamics of the subnetworks. Both A s (q), and B s (q) depend on the volume flow q, and can be represented in an affine decompositionWe approximate the input-output characteristics of each system s ∈ {0, .., c} by moment matching of the volume flow dependent transfer function H(q, ω) = C(ω1 − A(q)) −1 B(q) in frequency domain. Considering the volume flow q as a timevarying parameter to the advection of the energy density, local Galerkin projections V s e (q e ) are performed for representative volume flow vectors q e . These are picked using a greedy strategy and ultimately combined to a global projection V s using a singular value decomposition for each s ∈ {0, .., c}. Application of V s to (1-3) defines t...