Terrestrial biosphere models offer a comprehensive view of the global carbon cycle by integrating ecological processes across scales, yet they introduce significant uncertainties in climate and biogeochemical projections due to diverse process representations and parameter variations. For instance, different soil water limitation functions lead to wide productivity ranges across models. To address this, we propose the Differentiable Land Model (DifferLand), a novel hybrid machine learning approach replacing unknown water limitation functions in models with neural networks (NNs) to learn from data. Using automatic differentiation, we calibrated the embedded NN and the physical model parameters against daily observations of evapotranspiration, gross primary productivity, ecosystem respiration, and leaf area index across 16 FLUXNET sites. We evaluated six model configurations where NNs simulate increasingly complex soil water and photosynthesis interactions against test datasets to find the optimal structure-performance tradeoff. Our findings show that a simple hybrid model with a univariate NN effectively captures site-level water and carbon fluxes on a monthly timescale. Across a global aridity gradient, the magnitude of water stress limitation varies, but its functional form consistently converges to a piecewise linear relationship with saturation at high water levels. While models incorporating more interactions between soil water and meteorological drivers better fit observations at finer time scales, they risk overfitting and equifinality issues. Our study demonstrates that hybrid models have great potential in learning unknown parameterizations and testing ecological hypotheses. Nevertheless, careful structure-performance tradeoffs are warranted in light of observational constraints to translate the retrieved relationships into robust process understanding.